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ABSTRACT 

We investigate the population of z = 3 damped Lyman alpha systems (DLAs) in a recent 
series of high resolution galaxy formation simulations. The simulations are of interest because 
they form at z = some of the most realistic disk galaxies to date. No free parameters are 
available in our study: the simulation parameters have been fixed by physical and z = 
observational constraints, and thus our work provides a genuine consistency test. The precise 
role of DLAs in galaxy formation remains in debate, but they provide a number of strong 
constraints on the nature of our simulated bound systems at z = 3 because of their coupled 
information on neutral H I densities, kinematics, metallicity and estimates of star formation 
activity. 

Our results, without any parameter-tuning, closely match the observed incidence rate and 
column density distributions of DLAs. Our simulations are the first to reproduce the distri- 
bution of metallicities (with a median of Zdla — ^o/20) without invoking observation- 
ally unsupported mechanisms such as significant dust biasing. This is especially encouraging 
given that these simulations have previously been shown to have a realistic < z < 2 stellar 
mass-metallicity relation. Additionally, we see a strong positive correlation between sightline 
metallicity and low-ion velocity width, the normalization and slope of which comes close to 
matching recent observational results. However, we somewhat underestimate the number of 
observed high velocity width systems; the severity of this disagreement is comparable to other 
recent DLA-focused studies. 

DLAs in our simulations are predominantly associated with dark matter haloes with virial 
masses in the range 10^ < Afvir/M© < 10^^. We are able to probe DLAs at high resolution, 
irrespective of their masses, by using a range of simulations of differing volumes. The fully 
constrained feedback prescription in use causes the majority of DLA haloes to form stars 
at a very low rate, accounting for the low metallicities. It is also responsible for the mass- 
metallicity relation which appears essential for reproducing the velocity-metallicity correla- 
tion. By z = the majority of the z — 3 neutral gas forming the DLAs has been converted 
into stars, in agreement with rough physical expectations. 
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1 INTRODUCTION 

One of the most difficult and most important questions to ask of 
any cosmological simulation is whether, even if it matches some 
known properties of the observed Universe, the route via which it 
obtained those results is physically meaningful. It is tempting to ar- 
gue that, with the degree of parameter-tuning available to the mod- 
em simulator (stemming from our inability to maintain a sufficient 



dynamic range, uncertainty in gas physics and in particular star 
formation and feedback prescriptions), attempts to match a small 
number of observed properties can succeed without representing a 
physical route to that success. A sensible test of any suite of sim- 
ulations, therefore, is to scrutinize its predictions for observed re- 
lations which were not considered in the process of planning those 
simulations. Success or failure of the simulation to match such re- 
lations cannot be equated to success or failure of the simulation and 
its predictions as a whole, but it can lend weight in either direction. 

In this paper we will apply such an approach to the se- 
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ries of galaxy formati o n simulations most recen tly described in 
iGovemato et alj ( l2007h . ISrooks et al] ( l2007h and iGovemato et alj 
l l2008h (henceforth G07, B07 & G08 respectively). The z = out- 
puts of these simulations contain more realistic disc galaxies than 
have previously been achieved. In particular, the simulated galaxies 
form rotationally supported disks falling on the z = Q Tully-Fisher 
and baryonic Tully-Fisher relations, have a distribution of satellites 
compatible with local observations (G07) and have reasonable stel- 
lar mass-metallicity relations (B07). 

It should be noted, however, that in common with other galaxy 
formation simulations, the mass in the bulge comp onent of the G07 
galaxies is overestimated (see also lEke et al.l200lh . The problem is 
essentially one of angular momentum loss; it seems that to pr event 
this, both high numerical resolution jKaufmann et al.ir2007h and 
better models of feedback from supernova explosions are impor- 
tant elements (G08). The exaggerated bulges cause rotation curves 
to decline unrealistically over a few disk scale lengths to ~ 75% of 
their peak value. However, these problems appear to be shrinking in 
magnitude as resolution increases - and, regardless, the simulations 
under consideration form galaxies at z = which are as realistic 
as current computational and modeling power will allow, so that a 
detailed study of their properties during formation is justified. 

In this paper, we will investigate at z = 3 the predomi- 
nantly neutral gas which gives rise to damped Lyman alpha systems 
(DLAs). These are systems with column densities of H I in excess 
of 2 X 10^° cm~^, seen in absorption against more dist ant luminous 
sources (generally quasars); for a recent review see Wolfe et"al] 
( l2005h . The particular limit is historical, corresponding to the col- 
umn densities expe cted if the Milky Way were to be viewed face on 
( IWolfeetal.ll986l) . but simple physical arguments suggest it makes 
a convenient distinction between trace H I in the ionised intergalac- 
tic medium (IGM, below the limit) and clouds which are predom - 
inantly composed of Hi (above the limit, see IWolfe et alj|2005h . 
The latter clouds must absorb, in an outer layer, the majority of in- 
cident photons capable of ionising hydrogen {hv > 13.6 eV) and 
are therefore termed "self-shielding". 

The existence of an intergalactic ultraviolet (UV) field arising 
from the cumulative e ffect of external galaxies and quasars (e.g. 
iHaardt&Madaull 19961) plays many roles in understanding the state 
of these clouds - not only does it affect the ionisation levels in the 
optically thin transition regions, but it also contributes substantially 
to the heating budget via photo-ejection of electrons from atoms. 
Consequently, gas cannot cool to form neutral clouds without first 
collapsing in the presen ce of a grav i tational mass with virial veloc- 
ity of tens of kms"^ ( lReesl ll986'; Ouin n et alj|l996 ). suggesting 
such clouds are associated with dark matter haloes and hence pro- 
togalaxies. This rough physical argument is verified by previous 
simulations ( many of which are listed in Table [T}, although in the 
simulations of lRazoumov et all ( l2006h and Razoumov et al. ( 200^ 
DLAs often bridge multiple haloes, extending into the intervening 
IGM. (We discuss this issue further in Section |53] ) 

Irrespective of their physical nature, it is simple to show di- 
rectly that DLAs contain the majority of Hi over all redshifts 
2 > (e.g. I^tler 1987), which suggests they have an important 
role to play in the global star formation history. Curiously, how- 
ever, a wide range of diagnostics provide compelling evidence that 
the star formation rates in typical DLAs are small (< IMQyr"^). 
These include the low characteristic metallicities and dust deple- 
tions (for a review see Pettini 2006), the extreme rarity of de- 
tectable optical counterparts and the faintne ss of Lya emissio n (a 
large number of studies are summarised in IWolfe et al.ll2005l) . In 
DLAs with detectable molecular hydrogen (H2), the local UV can 



be estimated from pu mping into high energy rotational levels (e.g. 
ISrianand et al.ll200a) . results suggesting star formation rates com- 
parable to the present day Milky Way (~ 1 Mq yr~^). However, 
since few DLAs are associated with detectable H2 absorption, these 
measurements are not necessarily representative of the wider pop- 
ulation. The est imated cooling ra tes from the recently developed 
Cll] technique dWolfe et al.ll2003l) lead to star formation rate den- 
sity estimates of ~ lO~'^M0yr~'^ kpc~^, although the exact inter- 
pretation of these results is complicated by various assumptions in 
the method and the unknown area of a typical DLA system. 

A further constraint on the hosts of DLAs is given by 
kinematic information encoded in unsaturated metal absorption 
lines. Both the neutral gas (traced by low-ion transitions such as 
SillA1808) and any surrounding ionised gas (traced by high-ion 
transitions such as C IVA1548) may be probed. The exact rela- 
tionship between the hig h- and low-ion regions is not entirely cer- 
tain (e.g. lFox et al.ll20()7[ and references therein). Emphasis in our 
work, and most other studies, is placed on the low-ion profiles since 
these presumably reflect the kinematics of the gas giving rise to the 
DLA itself 

The earliest syste matical survey of DLA low- ion velocity pro- 
files was conducted bv lProchaska & Wolf d i ll 997h . who suggested 
that the observed kinematics could arise from a population of thick 
cold rotating disks with a distribution of rotational velocities sim- 
ilar to that observed in local disk galaxies. This is at odds with 
physical intuition and the prevailing Cold Dark Matter (CDM) hi- 
erarchical cosmogony in the sense that it requires the halo mass 
function, and hence power spectrum of fluctuations, to remain al- 
most unchanged over 10 Gyr between 2 = 3 and 2 = 0. A 
view more compatibl e with the standard model seems desirable. 
iHaehnelt et alj ( Il998t) were quick to apply numerical simulations 
of structure formation to show that a fiducial population of 2 = 3 
CDM haloes were not incapable of producing velocity profiles sim- 
ilar to those observed - however, some details have proved more 
problematic as we describe below. 

One way of quantifying the simplest property of the kinemat- 
ics is to assign to each DLA a "velocity width" 1190%, roughly mea- 
suring the Doppler broadening of any unsaturated low-ion tran- 
sition (see Section [T2I for a more precise definition). The veloc- 
ity width may be presumed to give an indication of the underly- 
ing virial velocity of the system responsible, although there is no 
guarantee that a particular sightline will sample the entire range 
of the velocity dispe rsion within the system; the simulations of 
IHaehnelt et"ZI i ll 9981) suggested that ^90% ~ 0.6uvir with a large 
scatter. 



Simulations such as those by Katz et al. ( Il996bl) ; 
iGardner et all ( Il997j, 1200 ll) ; iNagamine et al.l ( l2004ij) focused 
instead on the cross-sectional size of haloes as DLAs. Taken 
with a halo mass function (which throughout this work, we will 
produce for a ACDM "concordance" scenario) such trends can 
be used to predict the overall rate of occurrence of DLA systems 
- a prediction which, in most simulations, agrees roughly with 
observations. Including the results of Haehnelt et al. ( 1998), one 
can further attempt to predict the relative proportions of systems of 
differing velocity width. In general, simulations underestimate the 
incidence of high velocity widths; this can be traced to the majority 
of the cross-section being assembled from low mass (~ 10^ Mq) 
haloes. Similar difficulties are also encountered in varying degrees 
by semi-analytic models of DLAs. These are, of course, not 
suited to producing the details of the line widths but they can 
help to pin down which physical considerations affect them (see 
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Reference(s) 


Type 


SF 


lonization/RT 


Max Vol'i) 


Gas Res(2) 


js^dLZ cL di. \ lyyoD) 


SPH 


Mone 


Plane Correction'^' 


(on A/Tr^^^3 

yzz ivipc f 


1 r\S.2-\/\ 

iU ^VlQ 


Oaianer et al. ( lyy /a) 


SPH 


None 


Plane CoiTection'-'^' 


\AZ Mpcj 


1 n8. 2 1\ If 
lU M0 


oai oiiei ei ai. ( lyy / 












Haehnelt et al. (1998) 


SPH 


None 


Den. Cutt-l) 


n/a(5) 


IO^^'Mq 


Gardner et al. (2001) 


SPH 


Yes, weakFB(6) 


Plane Correction'^' 


(17Mpc)3 


10**'2Mq 


Cen et al. (2003) 


Eulerian 


Yes, withFBtf') 


Hybrid'^) 


(36 Mpc)3 


llkpc 


Nasainine et al. (2004a) 


SPH 


Multiphase/GW(**) 


Eq. Thin/MP(*) 


(34 Mpc)^ 


IO^'-'^Mq 


Naeamine et al. (2004b) 


Adpt Eulerian''-') 




Non-Eq. Live RT/post-processor''^''' 






Razoumov et al. (2006) 


None 


(8 Mpc)3 


0.1 kpc 


Nasamine et al. (2007) 


SPH 


Multiphase/GW(8) 


Eq. Thin/MP(**) 


(14 Mpc)3 


IO^'OMq 


Razouinov et al. (2008) 


Adpt Eulerian'^) 


Basic 


Non-Eq. Thin/post-processor' '^'^^ 


(45 Mpc)3 


0.09 kpc 


This work 


SPH 


Yes, withFBC''") 


Eq. Thin/RT post-processor'^^' 


(25 Mpc)3 


10-*'0Mq 



Table 1. Selected previous simulations of DLAs. Notes on methods: (^)The largest volume simulated for the study, in comoving units. '^^The best gas 
resolution achieved in the study, which may not have been achieved in the largest volume. For SPH (Lagrangian) simulations we give the smallest particle mass; 
for Eulerian simulations we give the finest grid resolution (in physical units at z = 3). ("^'UV background in optically thin hmit; sightlines post-processed 
using plane parallel radiative transfer and ionisation equilibrium. t^'UVB optically thin, but in post-processing all gas particles assumed fully neutral for 
number densities n > 10~^ cm""^. (^^This study used a re-sampling technique to study the high resolution dynamics of a limited number of haloes and did 
not collect cosmological statistics. FB = Feedback, i.e. the deposition of energy into the ISM due to supernova explosions. By "Weak" FB is meant thermal 
injection only, which is generally recognized to be insufficient (see Section fsTll (^'The optical depth for each cell was calculated and used to approximate a 
local attenuation to the UVB. '^^The multiphase method keeps track of the fraction of a gas particle in a cold cloud phase in pressure equilibrium with the 
warm medium; the cold clouds are assumed to be fully self shielded while the ambient medium is regarded as optically thin. Phenomenological galactic winds 
(GW) are added, causing bulk outflows from all haloes. (^'Adaptive Eulerian: the grids which keep track of gas properties are automatically refined as regions 
collapse on small scales. ''^"'A comprehensive treatment of radiative transfer was used by this paper, using eight angular elements in the live simulation and 
192 in a post-processing stage. See Sections I2land l3 . l1 



e.g.lKauffmannlll996l ;l Mailer et alJIlOOll ; Ijohansson & Efstathioij 
\2Gm . 

Overall, the extent of the implications of the velocity mis- 
match is somewhat unclear. Some authors have argued that a fun- 
damental difficulty with th e ACDM scenario has been uncovered 
jProchaska & Wolfell200lb . However, numerical modeling of the 
DLA population is intrinsically troublesome. The dynamic range 
of processes involved in making the final population is tremen- 
dous, with the relevant scales ranging from cosmological to stellar. 
Therefore we suggest that failure to match - or, indeed, success in 
matching - specific observations should be interpreted conserva- 
tively. We now describe some specific uncertainties in understand- 
ing the simulated DLA population, and outline how these are dealt 
with in our work. For comparison, Table[T]gives details of previous 
simulations and their approach to these problems. 



of parent halo mass (l Gardneretal]|l997j|bl) ; observations of low 
redshift DLAs appear to confirm such a view jZwaan et alj|2005h . 
Presumably one must allow a sufficient volume in a simulation to 
study a number of protocluster regions as well as maintain suffi- 
cient resolution to resolve low mass dwarf galaxies. Techniques 
involving extrapolation into unresolved regimes (such as fitting 
functio nal forms to the DLA cross-section as a function of halo 
mass: [Gardner et al.ll 19973 15) provide useful signposts but natu- 
ral WjTiust_be_j^gardedwith caution. Our simulations , and those 
bv lRazoumov et al] j2006h and Naga mine et al.l l l2004ah resolve all 
haloes of relevance only by separately simulating a number of 
boxes of varying size; thus a way of combining results from the 
various boxes must be found. In our case, this involves using the 
halo mass function to re-weight sightlines in calculating cosmolog- 
ical properties (Section [331 >. 



(i) Star Formation. Although 'Gardner et al. argue that 
star formation at z > 2 has little impact on the column densi- 
ties of individual systems, it is not a priori clear how the kinemat- 
ics and cross-section are affected by the supemov ae feedback. In- 
vesti gations in this direction have been made by iNagamine et af] 
( l200 4a b), using a phenomenological galactic wind prescription. In 
the present paper, our star formation is fully prescribed by physical 
models and z = observations, leaving no parametric freedom; 
see Section|2] 

(ii) Self-Shielding. DLAs are known to contain gas which 
is largely self-shielding (see above); thus the fiducial "uniform 
UV background" which is used in many cosmological simula- 
tions proves inadequate. We use a simple radiative transfer post- 
processor to coiTect for this (Section r3.ll ). We assess both the algo- 
rithm's reliability and the severity of neglecting radiative transfer 
in the live simulation in Section [54l 

(iii) Cosmological Sampling. The total DLA cross-section is 
thought to be evenly spread through many orders of magnitude 



The remainder of this paper is structured as follows. In Sections|2] 
and [3] we describe respectively the series of simulations in use for 
our study and our methods for extracting DLA sightlines and pro- 
ducing quantities representative of a cosmological ensemble. We 
give results for these quantities and their underlying relations in 
Section|4] Section|5]describes a number of consistency checks and 
runs with altered parameters which shed further light on the ori- 
gin of some of our relations. Finally, we summarise and discuss 
future work in Section [6] Two appendices contain technical detail 
for completeness. 

We adopt the following conventions. Except where speci- 
fied, all quoted measurements are given in physical units; where 
cosmological parameters enter calculations, these are based on 
the standard cosmology used for the simulations: Q,m = 0.30, 
fib = 0.044, = 0.70, (78 = 0.90, h = J/o/(100 km s^^) = 
0.70, Us — 1. We briefly investigated the effect of incorporat- 
ing the favoured parameters from the fifth year WMAP results 
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Tag 


(Mp.gas) 


(A/p,DM) 


e/kpc 


Usable Vol (comoving) 


Dwf 


1O*-°M0 


IO^-OMq 


0.15 


5 Mpc^ 


MW 


105.2 


106.2 


0.31 


50 Mpc^ 


Large 


106.1 


107.0 


0.53 


600 Mpc^ 


Cosmo 


106.7 


107.6 


1.00 


15625 Mpc^ 



Table 2. The simulations used in this work. The first column is the tag which 
we use to refer to each simulation. For all except "Cosmo", a subsample of 
the full box is simulated in high resolution; no results are taken from outside 
this region. The second and third columns refer respectively to the mean gas 
and dark matter particles within the region, the fourth to the gravitational 
softening length (in physical units) and the final column gives the comov- 
ing volume of the region. The separate boxes are generated from entirely 
different sets of initial conditions; the "Dwf" and "MW" simulations are 
designed to form respectively a dwarf and Milky Way type galaxy at 2 = 
while the "Large" and "Cosmo" boxes are more statistically representative. 

jPunklev et al.|[20o3) . but the overall differences are expected to 
be minor (Section[53]l. 



2 SIMULATIONS 

The simulations are successors to the galaxy formation simula- 
tions described in G07, with higher resolution and a more ph ysical 
star f ormation feedback prescription as described bv lStinson et alj 
( I2OO6) (S06, see also below). They a re computed using the SPH 
code Gasoline iWadslev et al.l2004h . and include gas cooling and 
the effects of a uniform ultraviolet (UV) background (following 
[Haardt & Madau 1996) in an optically thin approximation. We later 
perform post-processing to account for self-shielding effects (Sec- 
tion |3.1| ). We also ran test simulations which included an approxi- 
mate treatment of self-shielding within the live simulation (Section 
El. 

The SPH smoothing lengths are defined adaptively to use the 
32 nearest neighbours in all averaging calculations, except that a 
minimum smoothing length of 0.1 times the gravitational softening 
(e) is imposed. The star formation and feedback recipe (S06) is 
based on the algorithms originally described by Katzl jl992l) . In 
short, gas particles can only form stars if T < 30 000 K, rigas > 
0.1 cm~'^ and the local hydrodynamic flow is converging. The rate 
of star formation in such regions is assumed to follow the lSchmidj 
( Il959f) law (pstar oc Pgas) with index n — 3/2. This is a rather 
natural choice of n, since it implies the cold gas is turned into stars 
over some multiple 1/c* of the local dynamical timescale: 

^^=c../G^^. (1) 

Pgas at 

Following the evolu tion of each star particle consistently with 
the iKroupa etS] ( Il993l) IMF, a fixed fraction esN of the super- 
nova energy created at each timestep is deposited thermally into 
the surrounding gas particles. To emulate the physics of the pro- 
cesses responsible for distributing this energy to the gas, radiative 
cooling is disabled in particles within a radius which must be de- 
termined. Traditionally this can involve further parameters, but the 
S06 method obviates the need for these by modeling the physical 
processes according to a prescription based on blast wave mod- 
els (e.g. [chevalier 1974; Ostriker & McKee 1988): this sets the ra- 
dius of the local ISM affected as a function of the local density 
and temperature. The free parameters, consisting of the constant of 
proportionality in the Schmidt law (c*) and the efficiency of su- 



pernova energy deposition into th e ISM (esn) are tu ned such that 
isolated galaxy models match the lKennicuttI ( Il998bl) law and cold 
gas fractions observed in present day disk galaxies. This leaves no 
free parameters in the context o f this stu dy (setting c* = 0.05 and 
esN = 4 X 10^" ergs - see Stinso n et al.| 2006) but produces galax- 
ies which satisfy a wide range of observational constraints (see In- 
troduction). As part of the supernova algorithm, metals are also de- 
posited into the ISM (see Section 14.2.3 1 for more details), but note 
that our simulations do not yet contain any contribution to radiative 
cooling from metals (which is expected to be physically dominant 
below temperatures T ~ 10^ K). 

As in G07 & G08, the prescriptions are implemented in a fully 
cosmological context, based on a fiducial model with parameters 
given at the end of Section [T] The majority of our results are de- 
rived from four simulations (Table |2j. Three of these take advan- 
tage of the "volume renormalization" technique, in which progres- 
sively higher resolution is employed towards the central galaxy of 
interest (Katz & White 1993). The "Dwf" and "MW" simulations 
simulate the same region as the DWFl and MWl runs in G07 and 
G08, forming at 2; = a dwarf and Milky Way type galaxy re- 
spectively. The third simulation, "Large", is a compromise between 
high resolution and volume, while the fourth simulation, "Cosmo", 
is of a full (25 Mpc)^ box populated with gas at uniform resolu- 
tion. The smoothing length is constant in physical units for 2 < 8, 
and evolves comovingly for 2; > 8; its final fixed value for each 
simulation is given in Table |2] Otherwise, the physics in each of 
these runs is identical. 



3 PROCESSING PIPELINE 

The processing pipeline is based on SiMArfl an object-oriented 
C-l~l-/Python-based environment for analyzing simulations of arbi- 
trary file format with support for OpenGL real-time visualization. 
Such a framework allows us to understand our data rapidly, espe- 
cially since it includes support for stereoscopic glasses - giving a 
true 3D rendering of the data. It is written in such a way as to allow 
analysis of data from an interactive python session, hiding a large 
number of details from the scientific user. 

3.1 Self-Shielding 

We employ a basic radiative transfer post-processor to assess the 
effect of self-shielding. The simulation is divided into nested grids 
so that the lowest level cells contain close to 32 particles, to cor- 
respond with the SPH scheme of G07. We use a cosmic ultraviolet 
backg round (UVB) following a revised version of lHaardt & Madaul 
l[l996) (Haardt 2006, private communication); this is initially as- 
sumed to be present throughout the simulation volume. 

In each cell, using the radiation field described above, the ion- 
ization state of the hydrogen and helium gas is determined using 
an equilibrium solver algorithm based on lKatz et al.l ( Il99^ah . This 
requires the temperature and density of each SPH particle, derived 
directly from the simulation output. We keep the temperature, not 
the internal energy, fixed during this process. While not ideal, it 
is necessary to make some such assumption to determine the ther- 
mal state in post-processing. As a test we verified that keeping in- 
stead the internal energy fixed makes little noticeable difference to 
our results. On the other hand, shielding can significantly drop the 

^ www . ast . cam . ac . uk/ " app2 6/ siman/ 
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Deeper level ignored - 
not directly connected 
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Figure 1. A two-dimensional representation of our three-dimensional radia- 
tive transfer scheme, which involves building a nested grid and integrating 
optical depths from a cell along that grid to the edges of the box. A combi- 
nation of local fine spatial resolution and speed is achieved by traversing the 
grid at the lowest available level until this level is closed. At such a point, 
the algorithm jumps to a higher level in the grid. Subsequently, when lower 
levels become available they are ignored since these high density regions are 
not directly connected to the original high density region. This is primarily 
a computational simplification, but it also (in an admittedly ad hoc fashion) 
prevents long distance unphysical shadowing which would otherwise arise 
from our limited angular resolution - see text for details. 



heating rate, so that an assessment of the dynamical effect is essen- 
tial; we find it to be minor in terms of our final results, although 
uncertainties remain (see Sections [5.2l and [53t . 

With the ionisation state determined, the attenuation of the UV 
field due to the Hi, He I and He 11 ions is then calculated. This 
translates into an optical depth (as a function of frequency) for the 
cell. In subsequent iterations, the attenuation of the incoming UVB 
radiation is calculated by integrating the cell optical depths along 
the six directions defined by the orientation of the grid to the edge 
of the box. Because the grid is refined adaptively, this process is 
somewhat complicated by the need to "level-jump", i.e. move to 
higher levels of the nested structure when the lower levels run out. 
Note that the grid walk stays at the lowest level in any directly con- 
nected region of the origin cell (including when crossing bound- 
aries of higher levels). As it moves out of the high resolution re- 
gion, higher levels are selected as appropriate. The walk does not 
re-descend, even if lower levels again become available, instead us- 
ing the averaged properties of regions defined by the current level. 
This is primarily a matter of computational speed; however, it has 
the side-effect of preventing long distance shadowing which can 
arise in low angular-resolution radiative transfer mechanisms. Ad- 
mittedly we have not formulated this in terms of any limiting proce- 
dure yielding a well-defined physical calculation, but heuristically 
the averaging over larger solid angles as the integration proceeds 
away from the target cell is quite correct. We have illustrated this 
situation in Figure [T] 

This entire process is iterated over the full simulation, conver- 
gence being assessed by changes in the UVB field and ionization 



state between steps. We found that the system defined above had 
some oscillatory behaviour on its approach to convergence, which 
led us to introduce a damping term which in effect averages the 
optical depth between iterations. This results in faster convergence, 
but does not essentially alter the scheme described. 

Our self-shielding process is crude when c ompared to recen t 
radiative transfer codes (for a comparison see llliev et al. IHooi); 
however, we do not believe this makes our results obsolete: see 
Section [54] A further complication must be considered, which is 
the failure of the scheme as described to account for any dynamic 
effects of the changing ionisation fraction and heating rates. We 
have assessed the severity of this problem using a local attenuation 
approximation: see Section [J!4l 

Figure |2] shows a z = 3 map of the neutral hydrogen in a 
200 kpc (physical) cube centred on the major progenitor of a Milky 
Way type galaxy in our "MW" box. It is coloured such that DLAs 
appear in red and Lyman limit systems appear in green/yellow. The 
locations and virial radii (see Section U!2l below) of dark matter 
haloes exceeding 5 x lO^M© in virial mass are overplotted. 



3.2 Calculation of Halo and Sightline Properties 

We start by locating individual dark matter haloes in t he simula- 
tion u sing the grid-based code AHlQ ( Knebe et al. 200ll: [Gill et all 
We then mark a spherical region of radius rvir (the radius 
within which the density exceeds the mean density of the simula- 
tion at the given redshift by a factor 178, by analogy with spherical 
top-hat collapse models) as belonging to the halo of mass il/vir- We 
also record the masses in gas, stars and H I within the virial radius. 
The virial velocity is defined as usual, v^-^^ — GM^ir/rvir- 

In the resampled simulations, we immediately discard any 
haloes which have been contaminated by low resolution particles 
from the outer regions. We also discard any haloes with fewer than 
200 gas particles or 1 000 dark matter particles. We determined this 
limit empirically by examining the point at which halo properties 
diverged in simulations at different resolutions - for more details, 
see Section lSTI 

Line-of-sight properties can be calculated from particle-based 
simulations by projecting quantities onto a grid. For instance, pro- 
jecting all gas particles in a simulation onto the x - y plane al- 
lows the column density along the z direction to be estimated by 
summing the mass in a grid square and dividing by its area. This 
is the approach taken in previous SPH simulations of DLA proper- 
ties. However, in initial numerical experiments we found that sight- 
line properties in our simulations were not robust to changes in the 
somewhat arbitrary grid resolution. 

We have instead calculated all quantities using a true SPH ap- 
proximation; for details, see Appendix [A] This results in a vary- 
ing spatial resolution of sightlines which is automatically consis- 
tent with the simulation data. The minimum smoothing length al- 
lowed is 0.1 times the softening (given in Table |2}, meaning we 
can resolve spatial gradients over as little as ^ 20 pc in our highest 
resolution simulation, although a more typical effective resolution 
is ~ 200 pc. 

For our main results, sight-lines are projected through the sim- 
ulation in random orientations and at random sky-projected offsets 
from the centre of a halo up to its virial radius. We verified that 
extending this search area to twice the virial radius had no impact 



|http : / /www ■ aip ■ de/People/aknebe/AMIGAj 



© 2008 RAS, MNRAS 000, 000-000 



6 A. Pontzen et al. 




—200' ' ' ' ' ' ' ' ' ' — 'is 

-200 -150 -100 -50 50 100 150 200 

x/kpc 

Figure 2. The ^ = 3 neutral column density of H I in a 400 Itpc cube centred on the major progenitor to a z = Milky Way type galaxy (box MW). The 
colours are such that DLAs (logj^Q Nm/ cm~^ > 20.3) appear in dark red and Lyman limit systems (20.3 > logj^g Nm/ cm~^ > 17.2) appear in green 
and yellow. The circles indicate the projected positions and virial radii of all dark matter haloes with M > 5 X IO^Mq. All units are physical. A stereoscopic 
version and an animated version of this plot are available at www . ast . cam . ac . uk/ ~ app2 6/. 



on our results (this can also be seen directly from Figure |2l(- This 
confinement of our DLA cross sections is discussed in Section [575l 
For each sightline, we measure directly the column density in 
neutral hydrogen. If this column density exceeds the DLA threshold 
(A^Hi > lO^" "^ cm~^), it is added to our catalogue. If not, it is 
immediately discarded; but we keep track of the numbers of all 
sightlines taken so that we may calculate 

_ f J^dla\ 

0"DLA = CTsoarch i^j 

where ascarch ~ T^r^h is the search area, ntotai is the total number 
of random sightlines calculated and hdla is the number of such 
sightlines which exceed the DLA threshold. In this way, we obtain 
a representative DLA cross-section for each halo without assuming 
any particular projection. 

We also produce an absorption line profile for a low-ion transi- 
tion such as those of Si II. We assumed that the relative abundances 



of heavy elements were solar and that Si II was perfectly cou- 
pled to H I, so that for solar metallicity Mx /Mh = 0.0133 aii d 
n(Sill)/n(Hl) = n(Si)/n(H) = 3.47 x 10^^ jLoddersll2003h . 
In other words, given the metallicity Z of each gas particle, we 
take 7i(Sill) = 3.47 x 10""' (Z/Zq) n(Hl). Although an approx- 
imation, we found that the effect of relaxing the assumption of the 
Si II -Hi coupling was minor; see Section [?!2l 

Example profiles are shown in Figure |3] in which we have 
chosen four haloes and displayed five random sightlines from each. 
For the purposes of this plot, we choose one of SiIIA1808, 1526, 
1304 or 1260 according to which transition has maximum optical 
depth closest to unity. The plots are centred such that Aw = cor- 
responds to the motion of the centre of mass of the parent dark mat- 
ter halo. We also added, for display purposes, gaussian noise such 
that the signal-to-noise ratio (S/N) is 30/1. The range contributing 
to 1190% (see Section |4. 2. 2t is shown in Figure |3] by vertical lines 
at either end. Absorption arising in haloes with virial velocities 
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Af/kms ^ Aw/km s 



Figure 3. Example low-ion (Sill) absorption profiles from random sightlines thi'ough selected 2 = 3 haloes: in reading order, these have virial velocities 
of 694, 207, 92, 59 km s"'^ (making the first halo extremely rare in cosmological sightline samples, see Section |4~ll . The first is taken from our "Cosmo" 
box, the second is the "MW" major progenitor and the remaining two are smaller haloes from the "MW" simulation. Note that the velocity axes are scaled 
differently in the respective plots. The zero velocity offset corresponds to the motion of the centre of mass of the halo (unobservable in practice), illustrating 
a qualitative shift in behaviour from multiple clumps (higher virial velocities) to single clumps (lower virial velocities). The vertical dotted lines indicate 
the velocity offsets where the cumulative optical depth reaches 5% and 95% of its maximum value; ^90% is given by the difference in their position in 
velocity space. The pixel size is approximately 5 km s^^, although internally we use a higher 1 km resolution. For the purposes of this plot, the profiles 
are normalized to correspond to a definite transition, although this is not necessary for our computations. This chosen transition, along with the sightline 
Hi column density and metallicity, is indicated in each panel. Noise is added to simulate S:N = 30 : 1; again this is only for illustrative purposes and is not 
part of our pipeline. 
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Wvir ^ 150 km s^^ is often composed of multiple clumps whereas 
for smaller haloes there tends to be one main, central H I clump, 
moving with the centre of mass of the halo. Visual inspection of the 
different systems suggests that the former systems are less dynami- 
cally relaxed, presumably because they have formed more recently 
and have longer dynamical timescales; however, we did not verify 
this systematically. 



3.3 Cosmological Sampling 

We aim to make statements about the agreement or otherwise of our 
simulations with cosmological observations of DLAs and therefore 
need to construct a represen tative global sample o f absorbers. Our 
approach is related to that of lGardner et alj ( Il997ah . in that we cor- 
rect our limited sample by reweighting in accordance with the halo 
mass function. However, we emphasize that in our case this is to 
make allowance for our li mited statistics and com bine results from 
separate boxes, whereas in lGardneretal] ( ll997iJ) this approach was 
used to extrapolate behaviour into an unresolved low-mass regime. 

Consider any measurable property of DLA absorbers, p. The 
aim of the process discussed below is to construct the distribution 
function, d?N/AXAp. (We will adopt the custom of using the ab- 
sorption distance X, where AX/Az = Ho{l + z)^ /H{z): popu- 
lations with constant physical cross sections and comoving num- 
ber densities maintain constant line densities AN /AX as they pas- 
sively evolve with redshift.) Given the monotonic invertible rela- 
tions p(A/vir) and (7DLA(A^vir), onc has 



(3) 



where auLA{M) is the cross section of a DLA of mass AI and 
f{M) is the halo mass functior^, so that the physical number den- 
sity of haloes with mass in the range M — + A/ + 5AI is f{AI)SAl, 
and Al/AX = c/Ho{l + z)"^. However, for an ensemble of haloes 
of the same mass, ctdla will be scattered - although it may be pos- 
sible to define a mean value (Jdla- Similarly, p will in fact be scat- 
tered around some suitably chosen average p{M) for a given halo 
mass; furthermore even for a single halo most properties p will vary 
depending on the particular sightline taken through the halo to the 
distant quasafl. 

The easiest approach is to make the replacements p — > 
p(A/vir), (J — > a"(Afvir) in equation l[3}. However, even if signifi- 
cant deviation of p from p is rare, one may be interested in varia- 
tions of A^N/AXAp over several orders of magnitude - these rare 
fluctuations can therefore contribute significantly to the "tail" of the 
observed values. 

Accordingly, the method by which we perform the reweight- 
ing is non-parametric. Since at first this can seem a little obscure, 
we offer two descriptions. Below, we have described the method 
in a heuristic manner. In Appendix IbI we have outlined how this 
method can be derived by discretizing a well defined integral, 
which allows for an exact interpretation of our results should one 
be necessary. 

First, we bin haloes (in logarithmic bins) by their virial mass. 
Within each bin i, ranging over virial mass Mi —> A/i+i, one may 



^ We adopt the numerically cali brated version of the ISheth & Torme j 

halo mass function given bv lReed etal]j2006h . 
^ We regard both of these effects as stochastic scatter, although presumably 
a complete theory would account for the exact value of p in terms of a 
sufficient number of parameters pertaining to the halo and sightline. 
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Figure 4. The DLA cross-section of haloes which meet our resolution cri- 
teria in the Dwf (plus symbols), MW (dots), Large (cross symbols) and 
Cosmo (tripod symbols) boxes plotted against their virial mass. There is 
a (resolution-independent) sharp cut-off at Mvir ~ IO^Mq below which 
the cross-section for DLA absorption is negligible. Haloes with no DLA 
cross-section are shown artificially at l og^ ^ °"dla/ kpc^ = — 1. The fit to 
the equivalent results for two models in lNagamine et alj j2004al) is given by 
the dotted and dash-dotted lines (their models P3 and Q5 respectively). The 
major progenitors to the "MW" (Milky Way like) and "Dwf" (Dwarf type) 
z = galaxies are indicated. 



calculate a mean DLA cross-section ai and a total physical number 
density Fi : 



f{M)AM. 



(4) 



Mi 



The contribution of DLA systems per unit physical length 
from the bin i is FiOi, and consequently one may calculate: 

diV dl ^ 

i 

To investigate the distribution of properties observed we have 
two approaches. The simplest is to extract a representative set by 
discarding selected sightlines. This is ideal for comparing corre- 
lations between sightline parameters (e.g. velocity & metallicity, 
Section l4.2.3t . where the observed data sets consist of only 0(100) 
separate sightlines. For each halo h, we need to select a number of 
sightlines proportional to Wh = o"?iJi(h)/"-i(h) (cf equation JBlOt ) 
where ah is the cross-section of the specific halo, i{h) represents 
the mass bin to which halo h belongs and rii is the total number 
of simulated haloes in the mass bin i. The actual sightlines chosen 
from each halo are determined by a pseudo-random but determin- 
istic approach, for stability of results. 

This method is not ideal, however, because it throws away po- 
tentially useful information. In particular, when constructing dis- 
tribution functions for a property p, we use an alternative method 
which uses the cosmological halo mass function to weight, rather 
than select, results. The set of all sightlines is binned by the prop- 
erty p, indexed by j. One then has 
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d'iVDLA 



dXdp 



E 



Pj 

Wh(k) 

Ap 



Num. DLA obs. through h with p in bin j 
Total DLA obs. through h 



(6) 



where the sum is over all sighthnes from any halo in any box, h{k) 
denotes the halo associated with the sightline k and Ap is the bin 
size. The constants are not hard to calculate, but obscure the tech- 
nique and are presented in full in Appendix IB] equation l IBSt . 

Of course, the above methods require that line-of-sight effects 
(such as the coalignment of multiple DLAs along a single sight- 
line) are not a dominant effect. We verified that extending all our 
sightlines through the entire boxes made no significant difference 
to any of our results. This is because the DLA cross-sections per 
halo are small so that even when strong clustering is taken into ac- 
count, the probability of a double-intersection along a line-of-sight 
is very low. (Note also that velocity widths, Section |4.2.2| are al- 
ways determined from unsaturated line profiles, so that trace metals 
in the IGM are too weak to change the measured widths.) 

One should also require that environmental variations of halo 
properties (in effect systematic variations with parameters other 
than halo mass) are unimportant. This requirement is harder to ver- 
ify, but we neither detect variations of DLA properties of differ- 
ing regions in our "Cosmo" box, nor any strong correlations with 
indirect measures of environment such as the halo spin parame- 
ters. More systematic work in the form of the GiMIC project, res- 
imulating carefully c hosen regions of the Millennium Simulation 
jSpringelet ai]|2005l) . also suggests that environmental considera- 
tions are largely controlled by the local variation of the finite vol- 
ume halo mass function (Grain et al., in prep). 



4 RESULTS 

4.1 Masses of Haloes Hosting DLAs 

Taking all haloes of sufficient resolution (Section lSTt . we plot their 
DLA cross-section against virial mass in Figure|4] Our DLAs have 
cross-sections of between 1 and lOOOkpc^. In agreement with 
rough physical expectations, the cross-section increases as a func- 
tion of mass. In their regions of overlap, the results agree between 
boxes: although the Cosmo box has a vastly larger volume than 
our other boxes, and hence exhibits more scatter by probing rarer 
haloes, we verified that by restricting the number of haloes at a 
given mass the distributions of haloes from differing simulations 
were in agreement in each mass bin. Two extreme outliers in the 
Cosmo box were individually inspected and found to be systems in 
the middle of major mergers. 

A notable feature is that the cross-section drops off extremely 
quickly for Mvir < IO^Mq (i.e. v^i^{z = 3) < 30 km s"^). This 
behaviour is accompanied by a sharp break in the mass of neutral 
hydrogen in such haloes (Figure|6ll, and we attributed it to th e UV 
field which prohibits cooling of gas in low er mass haloes (see lReesI 
1 1 98d lEf stathioull 1 9921 : [Ouinn et al.ll 1 996l and Section[531l. 

Our cross-sections are overall some what larger than previ- 
ous simulation works have suggested (e.g. iNagamine et alj|2004al : 
[Gardner etal] 1 1997bl) and furthermore are not compatible with a 
singl e power-law link to the h alo mass. The results of two models 
from INagamine et alj ( l2004al) are overplotted on Figures |4] and |5] 
for comparison (runs "P3" and "Q5" here refer to weak and strong 
galactic winds in the cited work). Our enhanced cross-sections for 
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Figu re 5. The d ata of Figure \4\ multiplied by the halo mass function 
from i Reed et alj (2006) to give a total line density for each represen- 
tative system. Haloes with zero cross-section are shown artificially at 
logio d^N/dXd logio M = -4. 
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Figure 6. The total mass of neutral hydrogen plotted against the virial mass 
for our haloes. The symbols are as described in the caption of Figure|4] 



Mvir ~ 10^" M0 are apparently a consequence of our particular 
feedback implementation - see Section lSTI 

From our cross-sections, we can calculate directly the line 
density of DLAs Zdla = dN/dX via equation ([5]l. For this, 
we find luLA ~ 0.070, in good agreement with the obser- 
vationally determined value Idla{z = 3) = 0.065 ± 0.005 
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from SPSS DRjfl usin g the method described for SDSS DR3 in 
IProchaska etal] lioOSh . In Figure |5] we have shown how haloes 
of different masses contribute to this total line density by plotting 
d^iV/dXdlogio = (InlO) (d;/dX) A//(M)crDLA for each 
halo; /dla is simply the integral under the curve defined by the lo- 
cus of these points. The major contributors to the total line density 
are haloes of masses IO^Mq < Mvir < IO^Mq. At lower and 
higher masses the contribution is cut off by the rapidly decreas- 
ing cross-sections or exponential roll-off in the halo mass function 
respectively. 

As expected from our previous discussion, our results have a 
peak at ~ 10^'^ M0 which contrasts with the flatter results from 
fiducial power-law cross-sections. Consequently, at first glance, it 
appears that the area under our locus of points must be larger than 
that under the N04 curves and hence a substantial disagreement in 
line density is inevitable; however the cut-off for N04 is at rather 
lower masses (Afvir ~ 10* M©), and this brings the total line den- 
sity in N04 considerably closer to the observed value. 

Plotting the total Hi mass against the virial mass (Figure 
^ and DLA size against the total H I mass of a halo (Figure [7]l 
gives an alternative view of our cross-sections. A striking fea- 
ture of the latter plot is a bifurcation, particularly notable in the 
"Cosmo" box but also traced by the "Large" box, in which haloes 
of a fixed Hi mass A/hi ~ 10^ Mq can have different cross- 
sections. The primary physical distinction between the upper and 
lower branches is in halo mass: the former trace Hl-rich haloes 
with A/vir < 1O^°-^M0 and the latter trace a population of Hl- 
poor haloes with Afvir > 10^"'^ M©. Th is is reminiscent o f recent 
claims of bimodality in observed DLAs jWolfe et allboOSh . How- 
ever we limit ourselves, for the moment, to more general considera- 
tions, noting that the high mass, low a branch DLAs are extremely 
rare in our simulations (making up less than 2% of the total cross- 
section). Our method for generating cosmological samples (Section 
|3.3l l will propagate through any such bimodalities into our final re- 
sults without difficulty, assuming the Cosmo box has a representa- 
tive selection of haloes. See Section [63] for a further discussion of 
bimodality. 

A guide (Tdla - Mm relationship can be estimated by fixing 
a particular star formation rate. As explained in Section |2l in our 
simulations the instantaneous star formation rate is given by the 
Schmidt law. From this may be estimated a timescale for conver- 
sion of neutral gas into stars. 



With the assumption that pgas — A'/hi/(0.74 cy^^j^, one has the 
approximate relation 



O-DLA ~ 10 kpc' 




Our locus lies roughly along the line r — lO"'^ years; this 
is not unexpected, especially if our simulations are to match ob- 
servations that Q,i{z = 0) ~ S1dla(2 ~ 3)0 This suggests that 
a large fraction of the DLA cross-section should be converted to 
stars by 2; = 0, giving an upper limit of r ;^ O(lO^'^yr) to the 
star formation timescale r (unless star formation proceeds in rapid 

^ www . ucolick . org/ "xavier/SDSSDLA/DRS/ 
^ In other words, the total mass of stars at z = in a given comoving 
volume is roughly equal to the total mass of neutral hydrogen in that same 
volume at 2: = 3; see equation j9) and the ensuing discussion. 
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Figure 7. As Figure |4] except the cross-sections ai'e now plotted against 
the total mass of neutral hydrogen in each halo. Haloes with no DLA cross- 
section are shown artificially at logio cdla/ kpc^ = — 1. The dotted lines 
are of timescales for H I depletion through star formation; from top to 
bottom T = lO^*^", lO"'^ and 10^" yr These illustrate constraints on our 
locus of points by considering both (2 = 0) and the lifetime of a typical 
DLA (see text for details). 

discrete bursts, which is not the case in our simulations). Assuming 
DLAs are not short-lived objects, or achieved by very fine balanc- 
ing of rapid gas cooling and star formation, one would also expect 
r ^ 0{1/H{z = 3)) ~ 0(10^' Vr)- The constraint r > lO^r is 
obeyed, suggesting that this stable model is reasonable. 

4.2 Column Densities, Velocity Widths and Metallicities 

4.2.1 Column Density Distribution 

One of the best constrained quantities, observationally, is the neu- 
tral hydrogen column density distribution f{Nm,X). This is de- 
fined such that f{Nm, X)dNHi<iX gives the number of absorbers 
with column densities in the range A'^hi ^ A^hi + dA^'ni and ab- 
sorption distance X X + dX. Applying the reweighting method 
described in Section [33] to our sample yields an estimate for the 
cosmological column density distribution, shown by the solid line 
in Figure [8] This can be compared directly to the observed distri- 
bution given by the points with error bars, which are derived from 
SDSS DR5 (see previous section for an explanation). The matching 
of the normalization and approximate slope of the observed column 
density distribution can be seen as a genuine success of the simula- 
tions: we emphasize that no fine tuning has been applied to achieve 
this result. Furthermore, our results appear to have converged at the 
resolution of the simulations used (see Section lSTt . 

From the column density distribution, one may express the 
total neutral gas mass in DLAs in terms of the fiducial definition 

noLA{z)= / ° / fiNm,X)NmdNm (9) 

C/HIpcO J1020.3 (.m-2 

where lUp is the proton mass, pcfi is the critical density today, 
(1 — /hi) — 0.24 gives the fraction of the gas in elements heav- 
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Figure 8. The simulations' DLA column density distribution (solid line) 
compared to the observed values from SDSS DR5 (points with error 
bars, based on Prochaska et al. 2005; see main text for explanation). The 
dashed, dashed-dotted and dotted lines show the contribution from Afvir < 
1OS-5M0, 1O9-5M0 < Mvir < 1O" °M0 and A/vir > IO^-'^Mq 
haloes respectively (these are not directly observable distributions, but give 
guidance as to how our cross-section is composed). 



ier than hydrogen and A'max is an upper limit for the integration, 
which is discussed in the next paragraph. nDLA(z) gives the frac- 
tion of the redshift zero critical density provided by the comoving 
density of DLA associated gas measured at redshift z. (This is dif- 
ferent from the more natural definition of time-dependent ils which 
express a density at any given redshift in terms of the critical den- 
sity at that redshift. Only in the Einstein-deSitter universe will these 
definitions coincide.) 

Although the calculation should take A'^max = oo, this is 
not possible for the observational sample owing to the rapidly 
decreasing number of systems at the high column density limit. 
|Prochaskaetai] ( l200 5') discussed how different assumptions for the 
functional form of the column density distribution can lead to dif- 
ferent values of fioLA- The discrepancies are small for the two best 
functional fits to the observational data (a double power law or a 
Schechter function with exponential roll-off at high column densi- 
ties). However, these extrapolations are actually only constrained 
by a few points at high column densities; a more robust approach 
- albeit less physically transparent - is to calculate JIdla directly 
from summing the total neutral hydrogen in the observed sample 
of DLAs, which for the SDSS DR5 sample is roughly equivalent to 
using the upper limit A'^max ~ 10^^'^^ cm~^. 

Using this limit, we obtain f^DLA.sim ~ 1-0 x 10""^, which 
can be compared with the result from SDSS DR5 in the combined 
bin 2.8 <z < 3.5, r^DLA.obs = (0.84±0.06) x 10"^ . As expected 
from FigurejS] the results are in fair agreement; the slight mismatch 
is driven by the overestimation of our high column density points 
(lO^i-S < iVHi/cm-2 < 1021-75). 

Our weighting approach predicts the form of the distribution 
for much rarer, higher column density systems than the sample 
limited observations allow. Significant contributions to J^dla are 
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Figure 9. The simulations' DLA velocity distribution (solid line) compared 
to the observed values based on the sample described by Prochaska et al] 
12003) (see text for details; shown by points with error bars). The dashed, 
dashed-dotted and dotted lines are as described in the caption of Figure[8] 

made by these rare systems unless 7 = dln/(A'^, X)/d\'aN <^ 
—2. In fact, directly measuring the slope 7 for our simulations 
shows that it slowly decreases from 7 ~ —1.0 for A'^hi ~ 
lO^"'^ cm~2 to a constant value of 7 ~ -2.5 for TVhi ^ 
10^^'^ cm~2. Thus a correction to f2DLA,sim is expected if we al- 
low A'^niax to extend to arbitrarily high values. Performing the cal- 
culation with A'^max = 00 gives ilDLA.sim = 1-4 X 10~^. This 
value is not directly comparable to observational estimates, but 
shows that an observer living in our simulations would underes- 
timate Odla by about 30% due to missing contributions from the 
rare high column density systems. A further discussion of this issue 
is given in Section|6j4] 

4.2.2 Velocity Width Distribution 

We have already discussed some qualitative features of the low- 
ion velocity profiles generated in our simulations (Section [T2l Fig- 
ure [3]l. These are important because they provide a direct observa- 
tional measure of the kinematics of the DLAs, and therefore have 
the ability to substantially constrain the nature of the host haloes. 
We now turn to the comparison of our characteristic velocities with 
the observed quantitative distribution. 

We assign a velocit y width to each generated p rofile using the 
fiducial 1190% technique jProchaska & Wolfelll997h . This inspects 
the "integrated optical depth" r(A) = dA'r(A') and assigns the 
velocity width fgo% = c(A(, — Aa)/Ao where T{Xb) = 0.95r(oo) 
and T(Aa) — 0.05T((X)). The result is a representative velocity 
width for the sightline, produced without any of the difficulties as- 
sociated with fitting multiple Voigt profiles. 

To a good approximation, the only dependence on the partic- 
ular low-ion transition chosen is in the overall normalization of the 
optical depths from the relative abundances and oscillator strengths. 
The ^90% measure of the velocity width is invariant under such 
rescalings, so that only for display purposes (i.e. Figure [3j do we 
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Figure 10. The data of Figure |9] (thick sohd Une) plotted cumulatively for 
comparison w ith similar plots in the l iterature. The dotted and dash-dotted 
lines show the lRazoumov et al] j2008l) models NI and HI respectively (see 
text for details). The shaded band gives the approximate Poisson errors on 
the observational data, but we caution these are in fact strongly correlated. 
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Figure 11. The simulated distribution of metallicities (thick solid line). The 
points with error bars show th e observed me tallicities from the 2 < z < A 
sample of DLAs based on Prochas ka et alj (2007), normahzed to the ob- 
served ^DLA- The dashed, dashed-dotted and dotted lines are as described 
in the caption of Figure[8] 



need to choose a particular ion and transition. Observers require to 
choose unsaturated lines because, while our simulations calculate 
T directly, spectra only determine e^^ which cannot be inverted (in 
the presence of noise) for r ^ 1. 

Using our weighting technique, we calculate the cosmolog- 
ical distribution of velocity widths, (i? N /AXdv. This is plot- 
ted in Figure |9] (solid line), along with observational constraints 
(points with error bars) calculated from data provided by Jason 
Prochaska (private communication, 2008) based on the compila- 
ti on of high resolut ion ESI, HIRES and UVE^ spectra described 
in IProchaska et al] (12003). The dataset consists of 113 observed 
DLAs with 4.5 > ^dla > 1-6 (we directly verified that the ve- 
locity width redshift evolution over this range is negligible, so that 
we use the full range of observed systems to bolster our statistics). 
We normalize the line density of DLAs in the observational sample 
to match the observed /dla(z = 3) = 0.065 (Section r4.2.1| ). 

Overall, the simulations reproduce the approximate pattern 
of observed velocity widths, with a peak in the distribution at 
V ~ 50kms^^; the agreement is fair for velocity widths v < 
100 km s^^. However, it is now a familiar feature of ACDM simu- 
lations that they produce too few high velocity width systems (see 
Introduction) and our results are no exception. We provide a dis- 
cussion of this discrepancy in Section [63] 

Recent velocity w idth results by Ra zoumov et al] ( l2008h and 
iRazoumov et al.l jlOOd) have been presented by plotting a cumu- 
lative line density, Zdla(> v). A direct comparison between our 
results and those from such work can be drawn from Figure (TO] in 
which we similarly replot our velocity width distribution cumula- 
tively (thick solid line). As well as the observational data (thin solid 

^ Echellette Spectrograph and Imager on the Keck Telescope; High Reso- 
lution Echelle Spectrometer on the Keck Telescope; Ultraviolet and Visual 
Echelle Spectrograph on the Very Large Telescope 



lin e) we have overp lotted results from two simulations described 
bv lRazoumovetal.1 (2007; R08). These two models, R08.NI and 
R08.H1, differ in their size and hence resolution; the former covers 
a volume of approximately (5.7 Mpc)'^ (comoving), resolving grid 
elements of minimum side length 90 pc (physical) while the corre- 
sponding values for the latter are (46 Mpc)"^ and 0.7 kpc. Compar- 
ing our simulations with those of R08, the extent of the disagree- 
ment between simulated and observed results is similar (discount- 
ing the lack of low velocity systems in R08.H1, which arises from 
the coarser resolution), although our own results are actually some- 
what closer for low velocity widths {v < 100 km s~^). This is pos- 
sibly linked to our somewhat higher than normal cross-sections for 
haloes of virial velocities Uvir ~ 100 km s^^, but because much of 
R08's DLA cross-section lies outside haloes it is hard to provide a 
concrete explanation for such disparities (see Section|53}. 

We caution that any cumulative measure has strongly corre- 
lated errors so that the plot can present a somewhat distorted picture 
of the discrepancy (one is really interested in its gradient in linear 
space). In our plot we have shown the Poisson errors on the ob- 
served data as a shaded band, but these only represent the diagonal 
part of the covariance. 

4.2.3 Metallicity 

The metallicity, i.e. the ratio by mass of elements heavier than 
helium to hydrogen, is an important diagnostic of observed DLA 
sightlines. Since metals are deposited in the interstellar medium 
(ISM) through supernova explosions, the metallicity of a region 
is determined by the interplay between the integrated star forma- 
tion rate and bulk motions of gas including galactic inflows and 
outflows. In general, one observes a positive correlation between 
the mass of a galaxy and its meta llicity both at z = (e.g. 
iTremonti et Zll2004Eee et al.ll2006l) and at higher redshifts (e.g. 
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ISavaglio et alj2 005':'Erb et al."2006^. This trend is also observed in 
our simulations tBrooks et aU2007i . henceforth B07). Debate over 
a precise account of the origin of the relation continues, but the 
analysis of B07 shows that, for our simulations, the predominant 
effect is reduced star formation in low mass haloes, with the dy- 
namics of outflows providing an important correction to the simple 
closed box model. 

There are many uncertainties in simulating metallicities, 
which reflect not only the mass formed in stars but also the dynam- 
ics of the gas, its accumulation from the IGM and mixing within the 
halo. B07 showed that, for our simulations, high resolution is re- 
quired (A^DM > 3500) to attain convergence for the star formation 
history (and consequently the metallicities jf| . We verified (using 
an independent analysis code) that our own results suffered simi- 
larly, but that our other results converge at lower resolutions (see 
Section [STT l. Whenever we deal with results concerning metallic- 
ities, we therefore have to discard a number of haloes (those with 
A^'dm < 3500) which are included elsewhere. The effect of this 
is to reduce the number of haloes in each mass bin, and therefore 
exacerbate worries that we do not have a fully representative set of 
haloes for all mass scales. Nonetheless, we do not find any evidence 
that this is causing systematic effects, as we verified that our other 
distributions are not significantly changed by this restriction. 

Each of our sightlines is assigned a metallicity as follows. 

The simulations keep track of two sets of elements, the a-capture 

and Fe-peak elements. The enrichment process follows the model 

I . 1 I ' 

of 'Raiteri et al. (1996), adopting y ields for Type la and Type II 

supern ovae from Thiele mann et al.l ([l_986) and Weaver & WooslevI 
( Il993h respectively. For our purposes the differences between the 
two metal groups are minor: observations s how corrections are typ- 
ically < 0.3 dex (e.g. iLedoux et and exploratory work 
suggested our simulations were similarly insensitive. We choose to 
ignore these differences for the sake of simplicity and use the total 
mass density in metals. For quantitativ e results, we as sume a solar 
metallicity fraction of 0.0133 by mass ( lLoddersll2003l) . 

A cosmological distribution of DLA metallicities is generated 
using our fiducial technique (Section [3.3t . The result is shown by 
the thick solid line in Figure [TT] and closely matches the observa- 
tional constraint s, displayed as points w ith error bars. These are 
derived from the lProchaska et al] ( l2003h sample described in Sec- 
tion 14.2.21 restricted to 2 < 2 < 4. We use metallicities based 
on Q-capture element s (primarily Si), updated using data from 
IProchaska eTal] ( l2007h . 

Our simulated results are in good agreement with the ob- 
served distribution. They exhibit a roughly gaussian distribution (as 
a function of log metallicity); the best fit parameters give a median 
metallicity of [A//J/]mcd.sim = —1.3 with a standard deviation 
(T = 0.45. A similar fit to the observed data gives [M/H]^^^^^ = 
— 1.4 and a — 0.54. Given the uncertainties in abundances and 
yields, the differences are extremely minor. This success is unusual: 
previous simulations (e.g. Cen et al. 2003; Nagamine et al. 2 004b) 
tended to substantially overproduce metals in DLA systems, and 
were only been able to approach the observed result if metal rich 
DLAs could be hidden using substantial dust biasing. However, 
such a scheme is unsupported by observational evidence, and our 

^ This requireme nt is weaker than the analagous result in recent work by 
iNaabetaljfeoOTi) who require 10^ particles for convergence; t his probably 
relates to the neglect of feedback effects in lNaab et al.l j2007l) - feedback 
in our model tends to prevent collapse to very high densities, and thus im- 
poses an effective star formation scale independent of the resolution; see 
also Section lsTT] 
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Figure 12. The relationship between metallicity and low-ion velocity width 
of individual sightlines in a cosmologically weighted sample from our sim- 
ulation (crosses, with linear least square bisector given by the solid line) 
and from the observational dataset described in the text (dots, with linear 
least squai'e bisector given by the dotted line). The qualitative agreement, 
given the known deficiencies in the simulations and the lack of fine tuning, 
is remarkable. In the simulations, the origin of the correlation is an overall 
mass-metallicity relation; see text for further details. 



simulations now suggest it is unnecessary; for a further discussion 
see Section lMl 

We also investigated the relationship between the mean metal- 
licity of the gas within the virial radius of a halo and the spread of 
DLA metallicities derived from SPH sightlines through that halo. 
We found that the sightlines on average displayed slightly higher 
metallicities than the halo mean (by <^ 0.5 dex) with a spread of 
~ 1 dex. We interpret the offset by noting that DLA sightlines 
sample the high density gas in which star formation is occurring 
- even if the rate is low (Section I4.3t , presumably the local star 
formation preferentially enriches these. Some haloes showed a de- 
tectable radius-metallicity gradient, but in general this was shallow 
compared with the overall scatter. 



4.2.4 Correlation between Metallicity and Velocity Widths 

Observationally, there is known to be a correlation between the 
low-ion velocity w idth and the metallicity of a DLA sightline. 
lLedouxetal] ( l2006l) presented a set of observations showing a pos- 
itive correlation between these parameters at the 6ct level and sug- 
gested that the relation could reflect a mass-metallicity correlation 
analogous to that seen in g alaxy surveys; the result was confirmed, 
using a separate sample, bv lProchaska et al.l ( l2008t) . 

Our matching of the metallicity relation and near-matching of 
the velocity width distribution gives us confidence to attempt to 
probe this relationship in our simulations. We employ the Monte- 
Carlo sample generator version of our halo mass function correc- 
tion code (Section[33} to produce a sample of 64 coupled velocity 
and metallicity measurements, matching the size of the observa- 
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tional comparison sample described in Section 14.2.21 restricted to 
2 < 2 < 4 as in Section l433l 

The simulation results are siiown as crosses in Figure [72l witli 
the hnear least-square bisector fi0 given by the solid line. The ob- 
servational sample is shown by dots (the errors on each observa- 
tion are small compared to the intrinsic scatter), with the linear 
least-square bisector fit shown as a dotted line. These fitted rela- 
tionships are parametrized as logjQ Aiisim = 2.5 + 0.58 [M/H] 
and log^o A^obs = 2.7 + 0.53 [M/H] respectively. 

Although the normalization of the relationship in our simula- 
tion is somewhat different from the observational sample, we em- 
phasize that the slope and overall trend, as well as the mean metal- 
licity of our sightlines, are correct and that qualitatively the results 
are in agreement. Given that the metallicity distribution (FigurefTTt 
is in close agreement with that observed, we suggest that the quanti- 
tative disagreement arises from the previously noted underestimate 
of the velocity widths in our simulation (i.e. one should interpret 
the discrepancy in the relationships in Figure[T2]as a horizontal, not 
vertical, displacement). It is also apparent visually that the obser- 
vational results show a larger scatter about the mean relationship 
than does our computed sample. We tentatively suggest that this 
points towards an ultimate resolution of the velocity width issue 
which involves "boosting" the width of certain sightlines, perhaps 
by local effects such as outflows, while keeping the contribution by 
halo mass roughly as outlined in Section [43] (although this is not 
the only conceivable interpretation: see Section [673t . 

The origin of our relation between velocity widths and metal- 
licity is the u nderlying mass-metallicity relation, as suggested by 
iLedoux et al.l (j2006). We verified this directly, but it can also be 
seen from Figures |9] and [TT] the dashed, dash-dotted and dotted 
lines in each case show the contribution from haloes with Mvir < 
10^'^Mq, 1Q^'=Mq < Mvir < 1O"''^M0 and 10^°-^Mq < Mvir 
respectively. Both the sightline velocities and the metallicities can 
be seen to be a strong function of halo mass, and it is this fact that 
leads to the final relationship. 

4.2.5 Other correlations are weak 

The correlation between column density and metallicity or velocity 
is weak in our simulations, in agreement with observations: contrast 
the strong dependencies of the velocity widths and metallicities on 
the underlying halo mass with the weak dependence of the column 
density (Figure[8j. Given the historical confusion over the evidence 
for correlation between the column density and metallicity, we have 
investigated this aspect of our observational and simulated samples 
in a separate note (Pontzen & Pettini, in preparation). 

4.3 Star formation rates 

In the Introduction, we noted that star formation rates (SFRs) in 
DLAs are typically low (<; IMeyr"^) and yet these systems con- 
tain most of the neutral hydrogen, a necessary intermediary in the 
star formation process. We also see this behaviour in our simula- 
tions, and now describe how it comes about. 

We associate a star formation rate with each halo by in- 
specting, at z = 3, the mass of star particles formed within 
the preceding 10** years. For very low star formation rates (< 

^ The linear least-square bisector method estimates the slope of a relation- 
ship between X and Y as the geomet ric mean of the slopes obtained by 
regression ofX{Y) and Y{X); see .Isobe et alj il990l) . 
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Figure 13. The star formation rates in our haloes, plotted against their virial 
mass. Also shown is the observationally determined characteristic LBG star 
formation rate from lReddv et alj (2003), ™d a very approximate range of 
star formation rates for DLAs derived using the C II] cooling rate technique 
I Wolfe et al. 2003) - for more details and a list of caveats, see main text. An 
empirical power law fit is shown. Haloes with no star particles are shown 
artificially at logjQ M* = —4.5. 
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logioM^/Meyr-i 

Figure 14. Using our non-parametiic weighting technique, the distribution 
of star formation rates in DLAs within our simulation is shown, with the 
same observed ranges indicated as in Figure [T3] There is qualitative agree- 
ment between our estimates for DLA star formation rates and the observa- 
tional constraints. 
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Afp, gas/10* yr~^ < lO~^M0yr~^), there may have been no star 
particles formed in such a period; in this case, we extend the aver- 
aging period by a variable amount up to 10^ years so that we can 
resolve star formation rates down to ~ lO~'*M0yr~^. 

The SFRs of our individual haloes are shown in Figure [131 
the horizontal dotted line shows the star formation rate of an L* 
LBG galaxy. This is derived from iReddv et alj l l2007h . wherein 
is stated MAB(170q A) = -20.8. We adopt the conversion ratio 
of iKennicutj il998 j) . but first correct the luminosity for dust at- 
tenuation by a fac tor of 4.5, which is the mean correction given 
jReddy et alj|2007l . section 8.5). We then divide the final result by 
1.6 in order to consistently use the Kroupa IMF (Section |2](, for 
which a larger number of mas sive, hot stars are formed relative to 
the Salpeter IMF assumed by Kennicutll ( Il998al) . This gives a final 
characteristic LBG star formation rate of 39 yr~^ . 

We also plot a shaded band indicating the approximate SFR s 
obtained for DLAs using the C II] technique jWolfe et alj|2003h . 
This is intended to serve as a guide only ~ there are a number of 
uncertainties in the position of this band, including the intrinsic 
complexity of the observations and the conversion from a star for- 
mation rate per unit area (we have assumed a DLA cross-section of 
1 < a/ kpc^ < 100, which is our approximate range for the most 
common DLA haloes - see FigureslHandO. 

By evaluating 



p. = / dA/vir,/(Mvir)M.(MviO 



(10) 



(using a binning technique to discretize the integral) we derive a 
global star formation rate at 2: — 3 of ~ 0.2M(T)yr~^ Mpc~^ . 
The UV dust-corrected and IR estimates in iReddv et all bOOTh 
place this value between 0.05 and O.2M0yr~^ Mpc~^, but these 
are again sensitive to the IMF, assuming a Salpeter form for their 
main results: for a Kroupa IMF one should again reduce the rate by 
~ 1.6. This caveat should be borne in mind, but overall our results 
are not unreasonable and a detailed understanding of remaining dis- 
crepancies is beyond the scope of the present work. 

More importantly for our purposes, the star formation rate 
roughly scales as oc Ml{^ (see fit in FigurefTSt. Given the low 
mass end of the halo mass function is an approximate power law, 
M~^^''^ , the overall dependence of the integrand of equation dlOt is 
rather shallow and a wide range of halo masses contribute to the 
global star formation rate. Thus, there is no inconsistency in the 



view that a typical (Mv 



10 "Mq) DLA has a low star forma- 



tion rate but that the DLA cross-section as a whole contributes sig- 
nificantly to the star formation, and is largely converted into stars 
by z = 0. 

We generate a cosmological DLA cross-section weighted 
sample of these star formation rates (Section[33}. Figure[T4lshows 
the resulting distribution, d'^A'^/(dXdlogj^Q _A/star). As expected 
by the range of halo masses contributing to the DLA cross-section 
(Figure|4l(, the star formation rate in most DLAs is much lower than 
that in visible LBGs. Qualitatively, the simulated rate of DLA star 
formation is consistent with or fractionally lower than the observed 
star formation rate estimates from the ClI* technique (see above). 
Since we have not included upper limits in our approximate obser- 
vational band, this is an acceptable result. 

Finally, we investigated the total stellar mass accrued in our 
DLAs (the integral of the star formation rate). This quantity can be 
estimated in surveys of LBG by fitting their spectral energy distri- 
bution deduced from multi-band photometry. In Figure[T5]we have 
plotted the distribution of our DLA stellar masses and compared it 
with the estimated range of observed LBG's stellar masses, using 
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Figure 15. The distribution of stellar masses of our DLAs. The vertical 
band and vertical dotted line show respectively the range of and median 
value for observed LBG stellar mass from IShaplev et all ^200 11) (adapted 
for consistency with our IMF). As expected, most DLAs have formed a 
relatively small population of stars, consistent with their low metallicities. 



the results o f lShaplevetai] ( l200lh adapted as described above for 
our Kroupa IMF. Most DLAs have formed a relatively small popu- 
lation of stars (the distribution ranges over lO'' < A/^/Mq < 10^^ 
but peaks strongly at Af+ ~ 1O'''^M0), consistent with their low 
metallicities. There is a strong dependency on the underlying halo 
mass, which we have shown by plotting the contribution from dif- 
ferent halo mass ranges. This arises not only because the baryonic 
mass rises roughly linearly with the virial mass, but also because 
the star formation in low mass haloes is substantially suppressed 
by feedback (see also Sections [4.2.3l and l5. It . 

The simulated z — 3 DLAs are clearly being consumed very 
slowly, and it will be a matter of considerable interest to understand 
the ultimate fate of the gas, especially given the apparent realism of 
the galaxies at 2: = 0. We intend to address this question fully in 
future work; for the present, we note that over 80% of the neutral 
H I in the major progenitor at z = 3 has formed stars in our Milky 
Way type galaxy (box MW) by z — 0. However, this accounts for 
only 10% of the total stars; 10% of the stars were in fact already 
formed at z = 3, with the remainder coming from accreted gas and 
stars in satellites. 



5 CONSISTENCY CHECKS AND EFFECTS OF 
PARAMETERS 

In this section, we will discuss the effect of changing some details 
of our simulations and processing, including the resolution and star 
formation feedback prescription parameters. Since the effects are 
rather minor overall, some readers may prefer to skip directly to 
our conclusions (Section|6]l. The extra simulations used for the dis- 
cussions below are summarised in Table[3] 



© 2008 RAS, MNRAS 000, 000-000 



16 A. Pontzen et al. 



Brief Tag 




(A^p,DM> 


e/kpc 


Comment 


MW 


lO^-^ Mq 


10'='-2Mq 


0.31 


As Tabled 


T R 


lU 




U.Do 


Lower Resolution 


.LR.ThF 


106.9 


106.7 


0.63 


Thermal Feedback 


.LR.SS 


106.9 


106.7 


0.63 


Self-Shielding 


Large 


106.1 


107.0 


0.53 


As Tabled 


.LR 


107.0 


107.9 


1.33 


Lower Resolution 


.LR.ThF 


107.0 


107.9 


1.33 


Thermal Feedback 


.SS 


106.1 


107.0 


0.53 


Self-Shielding 


Cosmo 


106.7 


107.6 


1.00 


As TableE] 


.SS 


106.7 


107.6 


1.00 


Self-Shielding 



Table 3. A summary of the comparison simulations used. 
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Figure 16. As Figure |4] except we now plot only two simulations (MW, 
Large; dots and crosses respectively) and variants thereon - see Table|3]for 
a description of these. 



5.1 Resolution and Star Formation Feedback 

In this section, we describe a further suite of simulations which 
probe the dependence of our overall conclusions on resolution and 
feedback strength. Previous studies of DLAs have been shown to be 
somewhat sensitive t o the resolution (e.g. Nagamine et al. 2004a; 
iRazoumov et al ]|2006) and the main difference between our simu- 
lations and previous work lies in the feedback implementation, so 
that both these considerations are worth exploring. 

We described in Section [J!2l our resolution criteria, demand- 
ing both that TVdm > 1000 and iVgas > 200. This is quite 
a conservative limit, based on rerunning two of our simulations 
("MW" and "Large") at lower resolution (for details see Table 
(Sj. For haloes with coarser resolution than this stated limit, we 
found that the DLAs started to exhibit scatter about the locus 
defined by their higher resolution counterpart, and a slight ten- 
dency to underestimate the DLA cross-section (while overestimat- 
ing individual column densities). These resolution effects are rel- 
atively mild compared to some previously reported effects (e.g. 



[Gardner et all l ll997ij) . iNagamine et"ai] ( l2004ah . | Razoumov et al.l 
(2006|)). We suggest that this is because our treatment of feedback 
(see Section|2j, which leads to efficient energy deposition from su- 
pernovae in dense areas, helps to impose an effective scale-length 
below which the particles composing the ISM do not collapse fur- 
ther. Ultimately, given the impossible task of maintaining sufficient 
dynamic resolution to track the actual gas components making up 
the ISM, this is not an unreasonable behaviour. Figure [76l shows 
the "MW" and "Large" simulations along with their low resolu- 
tion counterparts ("MW.LR", "Large.LR") and demonstrates that 
the simulations are in good agreement when restricted according 
to our resolution criteria. The velocity widths from these low res- 
olution simulations were also overall in agreement with their high 
resolution counterparts; however, as previously noted, the star for- 
mation histories and metallicities only converge for a more strin- 
gent resolution cut (Section|4j23)- 

We took the low resolution initial conditions and reran the 
simulations, replacing our normal feedback implementation with 
a purely thermal approach (i.e. the same energy per supernova was 
produced, but the cooling switch-off was not implemented). This 
has the well-known effect of causing the supernovae energy to be 
radiated away over much less than the dynamical timescale jKatj 
1992; Thacker & Couchman 2000) and is thus the weakest conceiv- 
able mechanism. It is almost certain to be unphysical, since the fast 
cooling times arise from the combination of high temperature and 
density. This is caused by averaging properties of unresolved re- 
gions; in fact the high temperature (blast interior) and high density 
(undisturbed exterior) regions are presumably separated in the true 
parsec-scale ISM. 

Our cross-sec tions arising from this approach lie considerably 
closer to those of INagamine et all |2004a) (Figure [76l). We also 
find more usual results for quantities such as the metallicity; be- 
cause our metallicities converge only at high resolutions (see Sec- 
tion 14.2. 3t quantitative comparisons between our low resolution 
runs are a little dangerous, but we see a very m uch weaker mass- 
metallicity relation (see also lBrooks et alfcoOTh , with metallicities 
all approximately ^q/5. This is, as expected, in poor agreement 
with observations but in better agreement with older simulations. 
Recall that our fiducial feedback formulation causes a reduction 
in metallicities primarily by suppressing star formation efficiently, 
rather than causing any bulk outflows (Brooks et al. 2007). 

We also note that the velocity widths arising from particular 
haloes in our thermal feedback simulations are comparable to, or 
somewhat (~ 5%) higher than, the velocity widths in our fiducial 
simulations. (The velocity widths agree between low resolution and 
standard runs.) This does not result in a closer matching of the cos- 
mological velocity width distribution, because the proportion of in- 
termediate and high mass haloes is reduced. We interpret this as 
suggesting that the origin of our velocity widths is chiefly gravi- 
tational, and that the thermal feedback runs form denser, compact 
H I regions. This is no surprise: there is no physics in our simula- 
tions that can plausibly give rise to large-scale bulk galactic winds. 
Future improvement in DLA velocity width modeling may there- 
fore arise from an understanding of what drives such flows, and 
whether cold gas can indeed exist within them. However this un- 
derstanding will augment, not replace, a sufficient feedback model 
for the local suppression of star formation; see Section [63l 

5.2 Si II and low-ion regions 

One possible explanation for the underestimated incidence of high 
velocity width absorbers in our simulations (Section 14.2. 2t is a 
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Figure 17. The probability distribution of the ratio r of the velocity width 
from the CLOUDY-based Si II ionisation runs to the original velocity widths 
which assume n(Si Il)/n(H l) = n(Si)/n(H). Although higher velocity 
widths do arise by considering these effects, the mean differences are not 
large enough to have a significant impact on, for example, the overall dis- 
tribution of observed DLA line widths (Figure|9). 



misidentification of the precise regions responsible for producing 
the low-ion profiles. Therefore, we briefly investigated the effect 
of relaxing our assumptions (Section [3.11 ) that Ti(Si Il)/7i(H l) = 
n(Si) /n(H). This assumption is based on comparing the energetics 
of the ion transitions; however there is no actual physical mecha- 
nism coupling the Si II stage to H I. This is to be contrasted with 
Ol, which is strongly coupled to H I v ia a charge transfer mecha- 
msm (e.g. lOsterbrock & Ferland|[20o3) . There are clear examples 
of DLAs in which the Ol width is significantly smaller than the 
Si ll width (see, for e xample, the lowest two panels in Figure 6 
of iLedoux et al.|[r998h and thus there must be contributions to the 
Si II velocity structure from outside the cold H I gas. 

Because DLAs are low metallicity environments (both obser- 
vationally and in our simulations), a first approximation to the sil- 
icon ionisation problem can be achieved by taking the radiation 
intensity from our hydrogen/helium radiative transfer code (Sec- 
tion [3TTJ and calculating the equilibrium state of Silicon. To ad- 
dress this in a simple way, we calculated a grid of Cloud'J"! mod- 
els which varied in the local density, temperature and incident ra- 
diation strength indexed by the single photoionisation parameter 
+7— .Hii+e. For each particle in our simulation, we calculated 
the value of Sill/Si(r, T, p) by interpolating these models. Then 
we recalculated our DLA sightlines (Section l3.2b using the new 
Sin values, rather than the old assumed values. We used exactly 
the same set of offsets and angles, so that the final reprocessed- 
catalogue can be compared on a per-sightline basis. Our sightlines 
extend through the entire box so that any trace gas outside our 
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haloes could theoretically contribute to the velocity width. How- 
ever, the combination of the decreasing neutral gas density and 
decreasing metallicities makes the intergalactic contribution to the 
unsaturated line widths extremely minor. 

We have plotted, for all sightlines, the probability distribution 
of the ratio r of the updated velocity widths to the original in Figure 
[771 A number of sightlines did, after reprocessing, display signifi- 
cantly higher velocity widths (as much as doubling in some cases). 
Thus for individual systems, these corrections can be important. 
However (noting that the j/-axis in the plot is logarithmic) it is clear 
that the significantly increased velocity width systems are too rare 
to make a significant impact on the distribution of velocity widths 
(Figure|9]l, a fact we verified explicitly. 



5.3 Cosmological Model 

The parameters used throughout this work were fixed when 
running the first of the simulations described. Since then, 
our knowledge of cosmological parameters has improved 
thanks to the expanding wealth of constraints from galaxy 
surveys measuring the baryon acoustic oscillations (BAO), 
supernovae surveys and most signific a ntly t hr ee- and five- 
year WMAP results ( ISpergel et ^ l200d iDunklev et"al] 
l2008l) . Our chosen parameters are (Om, fib, f^A, fg, /i, Ws) = 
(0.30,0.044,0.70,0.90,0.72,1.0), whereas the latest WMAP, 
BAO and SN results combined suggest (JIm, fib, fiA, o"8, h,ns) = 
(0.269,0.046,0.721,0.82,0.70,0.96). It is natural to question 
whether such differences can have a significant effect on our 
results. 

Given our confined cross-section (Section [4.1b , any effect can 
be split into two components: the change in the halo mass func- 
tion, and the change in the individual objects. Considering the halo 
mass function first, since at 2; = 3 the relative density of A to CDM 
is suppressed by a factor of 4^ — 64, the difference in fiA has 
a negligible effect. Further, the value of fi^/fiivi has only minor 
consequences for the transfer function (amounting to the form of 
the acoustic oscillations). The most important difference, therefore, 
should arise from the reduction of erg: this will decrease the scale on 
which fluctuations have become non-linear, and hence reduce the 
number of high mass objects. However, we also need to consider 
the lower value of . The high mass end of the halo mass function 
is less affected (since it is closer to the normalization scale), but at 
the low mass end (corresponding to high k) there are a lower than 
expected number of haloes. Overall then, the effect of updating our 
halo mass function is to slightly reduce the line density of DLAs (to 
^DLA = 0.049, which is marginally lower than the observed line 
density, but acceptable given the uncertainties in our simulations) 
and shift the cross-section to somewhat lower masses. However, we 
verified that this had negligible effect on, for example, our veloc- 
ity width distribution. This is because, fortuitously, the two effects 
produce an almost flat value of /wMAp(Af)//fiduciai(A#) ~ 0.7 
over 10* < M/Mq < lO". 

On halo scales, the ~ 5% shift in the value of fij, must be the 
dominant effect. We do not believe this will have a large effect on 
our simulations, being a small correction, especially since equilib- 
rium considerations determine the physics on these scales. A com- 
plete resolution of this question, however, awaits new simulations 
with the updated values. 
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5.4 Radiative Transfer 

We have incorporated a correction for self-shielding (Section 
13. It vyhich is coarse compared to some recent simulations 
( iRazoumov et alj200^ : lRazoumov et alj200l) . On the other hand, 
the complexity of the fine-structure of the ISM today (and presum- 
ably also at z = 3) means that even the most advanced simula- 
tions cannot cap t ure its " microscopic" (i.e. parsec scale) behaviour. 
iNagamine et al.l ( l2004al) employed a sub- resolution two-phase 
pressure equilibrium model f or the ISM dSpringel & HernquistI 
l2003l : lMcKee & OstrikeJ[T977l) . apparently obviating the need for 
radiative transfer since the cold clouds are assumed to be fully self- 
shielded and the warm ambient medium to be optically thin. This 
is a promising approach, but it is not entirely clear how it links to 
the messy observational picture in which atomic gas exists in dis- 
tinct warm and cold phases, ionized gas in warm and hot phases, 
and photoionized regions occur within the cold phase around star 
formation ( H II regions) as well as in the diffuse disk ISM. 

For this work, our essential aim was to identify large scale 
regions in which the cosmological UV background should be sig- 
nificantly suppressed. We verified our code using a variety of tests 
comprising comparisons of equilibrium ionisation front positions 
with simple constant temperature CLOUDY models. While our ra- 
diative transfer code performs well in simple plane-parallel setups, 
it resolves only 6 angular elements and its true 3D performance 
is therefore significantly degraded. However, we believe that for 
the purposes of the present study it is adequate: in both optically 
thin and optically thick limits, the angular resolution is of little im- 
portance; only in the transition regions will significant differences 
arise. We ran a final test in which we post-processed our "Large" 
box after rotating it through 45°, thus providing a different set of 
angular elements to the code. None of our DLA results were sig- 
nificantly affected by this change. 

We do not include the UV emission of star forming regions 
in our calculation. This would introduce a significant extra com- 
plexity to our algorithm without any real benefits for the following 
reason. The star formation rate of most DLAs is known observa- 
tionally - and in our simulations - to be < lM0yr~^ (see Section 
14.31 Figure[T4l. For such low star formation rates the local contribu- 
tion to the UVB is comparable to, or usually significantly smaller 
than, the 2: = 3 field. This can be estimated, for example, using 
the data of iBlackl 1987, or more usually is obtained directly from 
observational data - see Introduction. Although the theoretical pic- 
ture of the importance of local sourc es is not straightforward (e.g. 
lMiralda-Escudil2005l : ISchavell200d) . we are confident that these 
direct observational upper limits justify our approximation. 

There is another slight inconsistency in our approach, in that 
the radiative transfer post-processor identifies certain regions of gas 
as neutral which in the live simulations is treated as ionized; fur- 
ther, UV heating occurs for the uniform background even in regions 
which are later identified as self-shielding. To assess the magnitude 
of this effect, we re-simulated "MW.LR", "Large" and "Cosmo" 
using a local UVB attenuation approximation ("SS" simulations in 
Table |3j. In these simulations, for each gas particle the strength of 
the UVB used for ionization equilibrium and heating calculations 
was reduced by the mean attenuation for particles of that density in 
the post-processed simulation. We continue to apply our radiative 
transfer post-processor, but the corrections involved are smaller in 
the new outputs. 

In general, the differences between the "live" self-shielding 
and original runs were rather minor. They consisted of a slight in- 
crease in the DLA cross-section (~ 0.2 dex) and a reduction in the 
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Figure 18. The mean (sight-line weighted) metallicity of DLAs in haloes 
from the "Cosmo" & "Lai'ge" boxes (tripod and cross symbols respectively) 
contrasted with their live self-shielding counterparts (dots and plus sym- 
bols). The most significant effect of including live self-shielding, other than 
a slight overall increase in DLA cross-sections, is a decrease in star forma- 
tion rates in low mass haloes, causing them to display systematically lower 
metallicities. For clarity in this plot, we have only displayed the first 100 
haloes from the "Cosmo" boxes. 



star formation rates. Consequently the cold gas metallicity of these 
low mass haloes was seen to drop by up to 0.5 dex (Figure [Tsl l 
- this was apparently the biggest change caused. In fact, this could 
help to resolve the slight paucity of low metallicity ([M/H]< —2.0) 
DLAs in our simulations (Figure fTTT l. although without a fuller set 
of statistics we were not able to verify this directly. Importantly, 
however, we reproduced all our measured distributions replacing 
"Cosmo" with "Cosmo. SS" and "Large" with "Large. SS". The ef- 
fects of this were dominated by an increase in overall line density 
to /dla = 0.080 (see Section [4.2. U . This takes us further away 
from the observed result /dla = 0.065 ± 0.005. However, except- 
ing normalization, the distribution of properties was left almost un- 
changed. We believe these effects warrant further investigation in 
the future, but their direct effects seem to be fortuitously small. One 
should be aware that, as well as the self-shielding effect which re- 
duces the UV heating and hence equilibrium temperature and pres- 
sure in the ISM, other inaccuracies in the detail of the ISM (such as 
magnetic pressure) could easily contribute opposite effects of sim- 
ilar magnitudes; however their study is well beyond the scope of 
this work. 



5.5 Intergalactic DLAs? 

It is clear from Figure [2] that our DLA cross-section is confined 
to within the virial radii of our host haloes. This is in qualitative 
agreement with the majority of previous DLA simulations listed in 
Table[T] and is an assumption of all DLA semi-analytic models that 
we are aware of . How ever, it contrasts with the re cent results of 
iRazoumov et al.l j2006l) and lRazoumov et al.l j2008l) (R06 and R08 
respectively) in which as much as 50% of the DLA cross-section 
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resides outside the virial radii of haloes (see R08 Table 2 and the 
top right panel of R06 Figure 3). 

Although rough estimates suggest that DLAs should form in 
dark matter haloes, it is not unreasonable for DLAs to form in 
overdense surrounding regio ns. Given the requir ement for self- 
shielding, nni ^ 10"^ cm~^ jHaehnelt et all 19981) , one may com- 
pare the gravitational mass required to confine such gas (Mgrav ~ 
NrnkT/GmpUm) to the total gas mass enclosed in the DLA 
(~ mpTi^^ N^j), giving the ratio 

Ms^^2(^] ( VV 1 (11) 

Thus, for DLAs with column densities close to the lower limit, it 
may be possible for the self-gravity of gas, or slight dark matter 
overdensities, to obviate the need for a fully collapsed dark matter 
halo. This is especially true if (unlike in our simulations) gas is 
allowed to cool efficiently to temperatures <C lO^K. 

The bigger problem is in understanding how the gas cools 
to become neutral in the first place. One key difference between 
our simulations and those of R06, R08 is that the latter include an 
approximate algorithm to correct the temperature for shielding ef- 
fects, whereas we keep our temperatures constant during the radia- 
tive transfer processing (Section lJTt . While our approach is clearly 
an approximation, the latter may overcool gas since dynamical and 
feedback heating effects could easily become significant as the UV 
field drops. The only way to correctly resolve this problem is (as in 
R06, but not ROB) to include radiative transfer in the live simula- 
tions. 

It will be interesting to see how this issue and a physical un- 
derstanding of it develop with future generations of simulations and 
radiative transfer codes. But while it remains to be adequately re- 
solved, we should note that (a) the effect in R06/08 is resolution 
dependent, and therefore the physical status is not transparent; (h) 
our self-shielding runs suggested the effect was rather small (Sec- 
tion |5.4t . at least with the local approximation and (c) we have pro- 
duced a DLA cross-section which matches nearly all observational 
constraints, which lends some indirect reinforcement to the validity 
of our approximation. 



6 DISCUSSION AND CONCLUSIONS 
6.1 Overview 

We have investigated the oc currence of D LA s in a s eries of simu- 
lations (see Govemato et al. l2007l [20o3 and [Brooks et al. 2007 - 
in this text referred to as G07, G08 and B07) which produce galax- 
ies at 2: = with as near as currently possible to realistic physi- 
cal properties (see Introduction). These high resolution simulations 
include a physically motivated star formation feedback prescrip- 
tion with only one free parameter (the efficiency of energy depo- 
sition) which is set by observations of low redshift star formation. 
Thus, as well as providing information on DLAs themselves, this 
work is an independent cross-check of the formation process of the 
G07/B07/G08 galaxies. 

To produce cosmological statistics from a smorgasbord of 
boxes, we used a non-parametric weighting process, in effect cor- 
recting the halo mass function of the combined sample. This does 
not involve any fitting or parametrization, and hence no extrapola- 
tion: all the results presented in this paper are derived directly from 
resolved regions of simulations. 

The picture that emerges from our simulations is of a cos- 



mological DLA cross-section predominantly provided by inter- 
mediate mass haloes, 10^ < Mvir/M© < lO". These ranges 
are somewhat higher than many early simulations suggested (e.g. 
[Gardner et anil997bl. l200l'). and are close to the range of lO^ '^ <, 
Mvir/ M(T) 10^^ sugges ted by observations of DLA/LBG corre- 
lations jCooke et al.[[2006l) . Our DLAs form stars at low rates (pre- 
dominantly <, O.IM0 yr~^) and have metallicities in the range 
-2.5 < [M/H] < with a median of approximately Zq/20, 
both in agreement with observations. We also investigated the dis- 
tribution of total stellar masses of DLAs, finding them to be pre- 
dominantly spread over lO'* < M*/M0 < 10"', with a peak at 
~ IO^'^Mq. By 2: = the majority of the neutral gas forming the 
DLAs has been converted into stars, but during this time substan- 
tial merging complicates the direct identification of low redshift 
galaxies with high redshift absorbers. In plots showing overall halo 
properties we have marked the location of the major progenitors to 
our z — Milky Way like and dwarf-type galaxies (from boxes 
"MW" and "Dwf" respectively); they appear to have been fairly 
typical DLAs. 



6.2 Mass-Metallicity Relationships 

A tight relationship between dark matter halo mass and star forma- 
tion rate (FigurefTBt underlies a strong mass-metallicity relation; in 
Figure[72lwe have shown that this is manifested by a correlation be- 
tween DLA veloci ty widths and metallicities (for equivalent obser- 
vational results see [Ledoux et al.l2003 ; [Prochaska et al.l2008l) . Fur- 
ther, the same simulations produce a realistic mass-metallicity re- 
lation for < z < 2 galaxies (B07), suggesting that these simula- 
tions capture the metal enrichment of collapsed systems in a mean- 
ingful way. The DLA result is not significantly affected by gradi- 
ents within the disks of our forming galaxies, but the DLA sight- 
lines do preferentially sample regions of gas with higher metallic- 
ities (by factors ^ 2) than the mass-weighted halo mean. Further 
work on metallicities in these simulations is ongoing, in particular 
regarding the importance of metal diffusion - the metals are sim- 
ply advected with the SPH particles, which can lead to inaccuracies 
where spatial gradients are important. In connection with this work, 
it will be interesting to see whether the metal enrichment of the sim- 
ulations' less overdense regions also appears realistic, providing a 
connection to the importance (or otherwise) of outflows (Section 
[U. 

Recent work has suggested that the origin of the observed 
relationship between kinematics as measured by Mg II rest-frame 
equivalent widths {W^^'^^^) and metallicity of DLAs should not be 
relied upon as an indicator of an underlying variation with mass 
(Bouche 2008). However, this claim is not directly relevant to fidu- 
cial ^90% velocity width determinations which are measured from 
unsaturated or mildly saturated transitions; Mgll is a strongly sat- 
urated absorption line, and as such Wr^''^^ can easily be affected 
by trace gas in the outer regions of galaxies (or perhaps even the 
ambient IGM). 



6.3 Velocity Profiles 

Overall, we reproduce the spread of velocity widths seen in DLA 
systems (Figure [9]l and approximately match the distribution of 
low velocity width systems {v < 100 km s~^). However the long 
standing difficulty of underestimating the observed incidence rate 
of high velocity widths {v > 100 kms~^) persists in our sim- 
ulations - although the discrepancy is rather small compared with 
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older simulations and c omparable to that seen in the recent work by 
iRazoumov et alj ( l2008h . It is not entirely clear why simulations in 
general have encountered such persistent difficulty in this respect. 
There are essentially three possibilities: 

Firstly, it is possible that our assumption n(Si Il)/n(H l) = 
n(Si) /n(H) oversimplifies the identification of regions responsible 
for the low-ion profiles and leads to systematic underestimates of 
the velocity dispersion. We investigated this avenue briefly (Section 
|5.2| l, but the correction did not appear sufficient to resolve the dis- 
crepancy - although this could be sensitive to the description of the 
ISM on sub-resolution scales (see also below). We ran a brief test 
in which all gas was assumed to contribute to the velocity width 
- in this (completely artificial) scenario, the cosmological rate of 
high velocity widths is over-estimated, showing that the required 
motions are, at least, "available" in this sense. 

Secondly, perhaps the internal gas kinematics still require 
some correction. Although we did not, for this work, attempt a rig- 
orous decomposition, qu alitative inspection sugges ted that a mix- 
ture of disk kinematics i Prochaska & Wolfdfl^Sl) and merging 
clumps jHaehnelt et alj[l998l) are responsible for the final spread 
of velocities. We did not find that our feedback made any sig- 
nificant contribution via turbulent motions to the velocity widths 
(Section I5.lt - hardly a surprise, given that the coupling to the 
ISM is achieved entirely thermally. A better understanding of feed- 
back is likely to help the kinematical situation by inducing bulk 
outflows. It would be interesting to see how the simulations of 
iNagamine et al.l ( 12004 a') perform in reproducing DLA line profiles, 
since they include a phenomenological model of galactic winds. 
Promising progress in placing s uch winds on a more physi cal foot- 
ing was recently described bv ICeverino & KlvpinI | |2007|) . whose 
simulations reproduce such bulk motions starting from an explicit 
model for the ISM. 

Thirdly, it is possible that our cross-section is associated with, 
on average, haloes of insufficient mass; in other words, DLAs at 
the high mass end could be overly compact. This in turn could 
be connected to the more well-known angular momentum prob- 
lem in which simulated disk galax i es have un derestimated scale 
lengths jNavarro & Steinmetzll2000l ; lEke et ahlboOl) . Our simula- 
tions, as was noted in the Introduction, go some significant way 
towards a resolution of the issue for z = di sk galaxies by 
combining sufficient Jfeedback ( Stin son et a i] |2006h with high res- 
olution ( Kaufmann et alJi2007n : our elevated DLA cross-sections 
for intermediate-mass haloes are connected to this same feedback 
(Section |5.U . As the mechanisms preventing angular momentum 
loss are further understood, these same mechanisms may continue 
to increase the DLA cross-section for high mass haloes. On the 
other hand, our matching of the metallicity distribution (Figure [TT) 
provides an alternative joint constraint on the star formation history 
and mass of responsible haloes - and is in good agreement with ob- 
servations. Any change in the mass of haloes comprising our DLA 
cross-section would therefore need to be compensated by a change 
in our virial mass - star formation relation (Figure[T3). 

6.4 Column Density Distribution and Dust Biasing 

As we discussed in Section l4.2.1l we match the overall column den- 
sity distribution well, but somewhat overestimate the number of ob- 
served systems with A'^hi > lO^'^'^ cm~^. This causes us to slightly 
overpredict value of SIdla (we obtain f^DLA.sim = 1-0 x 10^"^, 
whereas the SDSS data yield ^DLA.obs = 0.8 x 10"^). The sim- 
ulated value of nDLA.sim is obtained by imposing an upper cut- 
off on the column density of systems contributing to the mea- 



sured Odla because stronger systems are too rare to be observed 
with current observational datasets. But in our simulations, these 
same systems actually contribute significantly to the cosmic den- 
sity of DLAs; when included in the census, they boost the total to 
f^DLA.sim = 1.4 X 10^''. Thus in our picture, "invisible" systems 
contribute about 30% to the Hi budget; but functional fitting of 
the observed data from SDSS suggests a much steeper cut-off at 
high column densities than our simulations, and hence that the ob- 
served value OpLA.ob B given above should be essentially correct 
jProchaska et alj|200S ). While the precise observational extrapo- 
lation is subjective and the fit is driven by a couple of bins with 
A'hi > 10^^'^ cm~^, there is a strong constraint in the lack of 
systems at high column density. For instance, only three systems 
with A^'hi > 10^^ * cm~^ are observed in the SDSS DR5 survey, 
whereas our simulation predicts ten over the same path length (a 
glance at the Poisson distribution function shows that this is a sig- 
nificant discrepancy). 

It is notable that host galaxies of gamma ray bursts (GRBs), 
which trace a sample weighted towards higher column densi- 
ties, routinely display column densities of A'^hi ~ 10^^ cm~^ 
( Jakobsson et al. 2006) - so that such systems certainly exist but 
have a smaller cross-section than our simulations suggest. It is pos- 
sible that we overpredict the number of high column density sys- 
tems because our simulations are inadequate in the responsible re- 
gions, which are cold and dense: one should form H2 but the com- 
plex p hysics responsible is not implemented in our code. Schay^ 
( l200lh suggested that this molecular cloud formation could deter- 
mine a characteristic cut-off for the column density distribution. 
On the other hand, observations suggest that DLAs harbour very 
little molecular gas, with fractions at the very most /hj ~ 10^^ 
l lLedouxetai]|2003h . Of course the typical DLA sightline could 
miss high density clumps within giant molecular cloud formations, 
which would mean the global fraction of H2 could be somewhat 
higher. It is therefore currently unclear how feasible this explana- 
tion could be. 

It has been claimed in the past (e.g. ICen et al" I l2003h that 
DLA dust obscuration effects should substantially reduce the num- 
ber of observed high A^'hi systems relative to their true abun- 
dance. However, this is hard to reconcile with the most recent 
observations; in particular, radio-selected samples (which are un- 
likely to be biased by dust considerations) show little if any 
evidence for inc reased high TV hi absor bers jEllison eTa i] |200ll; 
Ijorgenson et alj2 006; Ellison et alj2008h . Further, in a companion 
paper Jpontzen & Pettinii2008i) we show that there is only marginal 
statistical evidence for dust-induced obscuration in the joint distri- 
bution of observed DLA column densities and metallicities, which 
are in fact nearly uncorrelated. This also arises quite naturally in 
our simulations, because the dependence of TVhi on the underlying 
halo mass is weak (Figure [8j, whereas metallicities and kinemat- 
ics are strongly correlated with the halo mass. Combining all the 
above considerations with our successful reproduction of the ob- 
served distribution of metallicities (Section|4j23J, we believe dust- 
biasing is unlikely to have a major part to play in a future under- 
standing of DLAs. 

The very slight trend which is observed linking the column 
density to the halo mass in our simulations biases low column 
density observations in favour of finding low mass haloes (Fig- 
ure [8]l. We have not explicitly studied the sub-DLA population, 
so that these resul ts are not directly comparable with the work of 
iKhare et al.l ( |2007|) . in which it is suggested that sub-DLAs pref- 
erentially probe more massive haloes than DLAs. However, there 
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is certainly a tension between these results which deserves some 
attention in the future. 



6.5 Bimodality 

Although we did not specifically set out to inve stigate the possibil - 
ity of bimodality within the DLA population jWolfe et al]|2008h . 
we saw a bifurcation in our cross sections as a function of H I mass 
(Figure |7](. The bifurcating lower branch in Figure [7] has higher 
virial masses for a given H I mass, wider velocity profiles, higher 
star formation rates and higher mean temperatures. It is important 
to emphasize that these systems, which are perhaps undergoing 
a starbursting phase, account for only ~ 2% of the DLA cross- 
section; our results are therefore not directly comparable to the ob- 
servational claim of two populations of roughly equal importance. 

We ran some careful tests to try to understand the origin of this 
effect, but it essentially remains work in progress: further simula- 
tions and detailed studies of the evolution of these objects will be 
key to a full account. The critical halo mass, Mvir ~ 10^"'^ M©, 
is probab ly too small to be directly related to the shock-stability 
model of iDekel^ Bimboid 12006) . There are no spatial correla- 
tions which would suggest the abnormal haloes are associated with 
protocluster regions. We verified that the results from the "Large" 
box were consistent with being drawn from a subvolume of the 
"Cosmo" box. Since the low cross-section branch objects account 
for such a small fraction of our DLAs, we felt justified in postpon- 
ing their full study to later work. 

6.6 Future Work 

In general, we have commented throughout that the lack of res- 
olution on fine ISM scales (a common feature of all current cos- 
mological simulations, likely to persist into the foreseeable future) 
is a fundamental limitation. The ISM is a complex mixture of gas 
in different phases with variations on truly tiny scales. It is also 
interesting to note that observations suggest that the pressure con- 
tributed by magnetic effects in the disk of our galaxy and other 
local gala xies is comparable in mag nitude to turbulent kinetic pres- 
sure (e.g. iHeiles & Crutche3 l2005'): hence magnetic effects could 
contribute at least as much as stellar feedback to an understand- 
ing of the formation of disk galaxies, at least at jz = 0. It is not 
impossible to fathom that magnetic fields could have an important 
role to play in a better understanding of DLAs. But regardless of 
any particular physical processes, our feeling is that long term fu- 
ture focus should remain on improving our understanding of the 
coarse-grained behaviour of such a mixture, i.e. in developing sub- 
gri d physics model s . This conclusion was reached independently 
bv lRazoumov et alj j20o3) . 

Even with our current set of simulations, it will be interesting 
to investigate in more detail the predictions of these simulations 
with regard to the LBG population; in Section 1431 we showed that 
the total star formation rate at z = 3 may be very slightly overpre- 
dicted compared to observation. But, in work not presented here, 
we have also investigated the shape of the luminosity function (at 
z = 3 and higher redshifts) and found both and the faint end 
slope to be consistent with observations. Remaining shortcomings 
could easily be accounted for by details of the IMF, suggesting that 
the underlying populations are quite reasonable (Govemato et al, in 
prep). 

In terms of the properties of DLAs, this work opens the door 
to an array of further studies. In particular, we intend to address in 



more detail the time evolution of our DLA population and the trans- 
fer of gas through hot, warm and stellar phases in the near future. 
Given the realisti c nature of metals i n our z = ?> DLAs and in our 
galaxies at z = jBrooks et al.l2007h , the details of the intervening 
time could, a mongst other t hings, shed light on the "missing met- 
als" problem ( |Pettini|[200^ . Further, it is known observationally 
that the metallicities, velocity widths and column density distribu- 
tions evolve only slightly with decreasing redshift; will our sim- 
ulations reproduce such slow evolution, which depends on finely 
balancing the cooling of gas, forming H I, with the formation of 
stars, destroying it? If so, the manner in which this is accomplished 
will be of considerable interest. (For instance, do individual ob- 
jects maintain the same DLA cross-section, or is the slow evolution 
only seen over cosmological averages?) Although flawed in many 
respects, simulations offer a guide to our understanding of these 
issues which should not be ignored - as long as it is taken with a 
healthy pinch of salt. 
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APPENDIX A: SPH CALCULATIONS 

For convenience we will denote the SPH averaging procedure by 
square brackets, i.e. for any quantity Q 

[Q\{r) = j A^r'W{fy-h{r'))Q{r') (Al) 
~ y^9iw{r,fi-hi) (A2) 

where W is the smoothing kernel and the smoothing length for each 
particle is given by hi, which notionally corresponds to a smooth- 
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ing field h{r). For each sight-line, the z axis is aligned such that the 
line of sight is parametrized by r = (0, 0, 2). The column density 
is calculated as: 



Nhi 



W^{\n\r,hi) 



Az [nm]{r) 



(A3) 



■ W^{\ri\^;K) (A4) 



i,rd.m(^^ + ^W (A5) 



where nmir) is the number density of H I ions, and nHi.i, rrn, pi 
refer to each SPH particle's H I number density, mass and density 
respectively and we have assumed a fiducial form for the kernel 
W{fl,f"2;h) = W(|fi - r2\/h)/h^. The projected kernel l |A5l l is 
calculated numerically; however, since conventionally W has com- 
pact support, so too does W±_ which need only be tabulated once for 
a finite range of values. The upper limit depends on the choice of 
kernel, but is conventionally max(d) = 2h. The sum iA4\ is taken 
over all particles, although in practice many particles are culled by 
considering their distance from the line of sight. 

In general, by projecting the kernel in this way one obtains 
an efficient and justifiable method for calculating any projected 
quantity. Integrated quantities follow directly by replacing nm,i 
in equation l lA4t with the appropriate particle property, subject to 
the following caveat: there is an ambiguity in calculations of this 
type involving more than one particle property or functions of par- 
ticle properties. In particular we will consider the generation of an 
absorption profile: 



'rtot(A) 



dzn{r)r.{X[l~v4f)/c];T{r)) 



(A6) 



where rtot(A) gives the total line of sight optical depth for wave- 
length A and Tv gives the absorption (Voigt) cross-section with pa- 
rameters assumed for the ion under consideration and temperature 
broadening T{f). There are (at least) two different methods for ap- 
proximating this. 

The first method is easiest to evaluate, and lumps the integrand 
into a single quantity which is to be SPH-interpolated: 



'rtot(A) = 



dz [nT,{X{l-v,/c)-T)]{r) 



(A7) 



^ — ^ /I, 



The second is considerably more computationally intensive, and in- 
terpolates each fluid quantity within the integrand before perform- 
ing the integration: 



'rtot(A) = 



dz [n]{f)T. (A (1 - b.](r)/c) ; [T]{r)) (A9) 



In the continuum limit W ^ S both these expressions reduce to 
iA6\ : further, it is simple to write a number of different expres- 
sions midway between which have this same property. For quan- 
tities varying on scales S> h the two methods will give identical 
results; however it is not clear a priori which, if any, method is op- 
timal for the DLA problem where as much information from near 
the resolution limit needs to be extracted. 

This problem is not confined to line of sight integrals, and so 
we may turn for guidance to standard approaches to SPH. Although 
we have not seen the ambiguity discussed explicitly, in general ex- 
pressions are used where the interpolation is performed as a final 



step as in dASt . We therefore adopt this approach; however we note 
that it has some odd properties, and in particular one may define a 
dispersion in any quantity 



\r) = [v\r} ~ [v]\r) ^ 



(AlO) 



with equality holding if and only if u(r) is constant within the 
smoothing region. In other words, the standard method forces us 
to accept a subresolution velocity dispersion in the fluid. The ex- 
act interpretation of this is somewhat elusive. Given a strictly la- 
grangian interpretation, at any point there is a well defined velocity 
and it should only vary over resolved spatial scales. But owing to 
the method of interpolation, if we find nearby two SPH particles on 
top of each other with equal and opposite velocities, this would not 
lead to any significant spatial variation mv{f) ~ but would cre- 
ate a non-zero [v]'^. The evolution of this "turbulence" would make 
for an interesting study, not least to determine its physical meaning, 
but is beyond the scope of this paper. 

Our final expression for the line profile is therefore 



Ttot 



(A) = r 

J — ; 



dzn,a{\{l-[v4)/c;T„{[v^] - [vf)(r)) (All) 



where the indicated thermal and turbulent broadening components 
are added in quadrature. Because is not a projected quantity, 
the z integration cannot be performed by simple integration of the 
kernel, and therefore this quantity is performed as an explicit nu- 
merical integration over the sightline. 

Note that this intrinsic velocity dispersion should be distin- 
guished sharply from the line of sight velocity dispersion, e.g. 



/ \2 , 2s _ f Jdz[nHlVz] 



J dz[nmv'i] 
J dz[nm] 



(A 12) 



It is clear, given the fiducial positioning of the SPH averag- 
ing brackets [] that this quantity will include contributions from the 
intrinsic velocity dispersion, as well as genuine line of sight disper- 
sion effects. 



APPENDIX B: WEIGHTING METHOD AS 
DISCRETIZATION OF AN INTEGRAL 

Let us consider a calculation of the line density of DLAs for the 
situation described in Section \3l2\ in which a has some range of 
values with probability densities P{a\AI) for a given halo mass 
M. Our starting point is the integral 



dM / daaf{M)P{a\M) 



(Bl) 



where we have temporarily adopted the physical distance measure / 
to obviate the need for dl /dX on the right hand side. We discretize 
this equation by assigning all haloes to bins in virial mass. Indexing 
these mass bins by i, and denoting the set of all haloes within each 
mass bin by Hi, one makes the replacement J dM f{M) . . . ^ 
• ■ ■ where Fi is defined by eq. ©. 
Our probability distribution for a, P{a\M), is represented by 
individual haloes. The simplest way to represent P is therefore by 
a sum of 5 function.'^ ~ one for each halo - normalized such that 



One could produce a more sophisticated method in which the 5 func- 
tions are replaced by smooth functions with unit area to create a smoother 
final P((t|A/). However, it is not clear to what extent this is helpful and 
both methods reduce to the correct integral in the iniinite-halo limit. 
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the total probability / P{(j\M)A(j 
1 



p((t|mo 



diVo 



dl 



E 



herii 



■ E 



(B2) 



(B3) 



where n(- • •) counts the elements of its parameter set, and Oh de- 
notes the cross-section of the particular halo h. This expression re- 
duces, as expected intuitively, to 

dA^DLA 



d; 



(B4) 



where Oi is the mean DLA cross-section of a halo in bin i. 

We are next interested in computing, for some property p, 

d'iVDLA 



AlAp 



^ da dMan{M)P{aSzp\M). 



(B5) 



This will involve exactly the same replacements as before, with the 
only difference being 

S{ah — cr) ■(T-^ 5{pj - p) 



P{ahp\M) 



E 



■E 



n{Sh 



(B6) 



Here, Sh is the set of all DLA sightlines taken through halo h. 
The expression 1 IB6I 1 will leave us with a floating S function in the 
chosen property p, so the final stage is to bin by this property. In 
other words, one integrates both sides of the equation with respect 
to p over a bin centred on pj with width Ap: 

^^' = E;^E-- 

Pi i ^ ' heHi 



dldp 



(B7) 



in which we have introduced Vj.h, the set of all sightlines in halo 
h such that p falls within property bin j (pj — Ap/2 < p < 
Pj + Ap/2). Note there is no sum over sightlines in the above (only 
a counting of sightlines in particular sets); the two sums are respec- 
tively over the mass bins and the haloes within them. Since every 
halo falls into exactly one mass bin, this is really just a sum over all 
haloes: 

j2. 



d'iVDLA 



dldp 



E 



Pi 



n{Hi(h)) n{Sh)Ap' 



(B8) 



where i{h) denotes the mass bin i to which halo h belongs. 

For completeness, with all constants (including the conversion 
to absorption distance X), our final calculation reads: 



sightline. From equation iB3\ . it is clear that the contribution of 
halo h to the line-of-sight density of systems is given by a proba- 
bility Wh (one may verify that this remains true for correlations be- 
tween any number of sightline properties) and thus the probability 
of accepting that sightline should be proportional to this quantity. 



d'iVoLA 



dXdp 



E 



Wh n{Vj,h) . 



Pi 



Ho{l + zr ^Ap n{Sh) 



F,, 



Wh 



(B9) 



(BIO) 



n{Hi(h)) 

There is no essential barrier to investigating correla- 
tions between properties by extending this method to produce 
d^A^DLA/dXdpdg for some extra property q. However, the grow- 
ing sparsity of the data in the increasing dimensional binning is 
such that an easier method suggests itself. Correlations between 
parameters in DLA studies are generally known only for 0(100) 
observed systems. One can generate from our simulations, using 
a Monte-Carlo technique, just such a limited sample; this can be 
compared directly to observational data. The probability for accept- 
ing a particular input sightline through this technique must be pro- 
portional to the output line density "generated" by that individual 
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